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ABSTRACT 


A novel numerical technique of solving convection-diffusion 
type equations, called the operator splitting technique has been 
used in the present study to solve model one and two dimensional 
heat transfer problems. Results of this study have been compared 
to those of two other conventional differencing schemes namely 
the upwind technique and the central difference technique. These 
results are further compared with the analytical solution 
wherever possible. It is found that for low Peclet number flows 
all the three schemes igive results identical to the analytical 
solutions. But at high Peclet numbers, convection dominates over 
diffusion and upwinding shows significant errors. Central 
differencing is also found to be unsatisfactory. However, the 
operator splitting solution comes close to the analytical 
solution. Hence we conclude that this scheme is relatively free 
from the maladies afflicting the conventional lower order finite 
difference schemes for a wide range of Peclet numbers. 

In the second part of the thesis, two phase oil— water flow 
through the porous media is numerically studied. In this 
problem, the operator splitting algorithm is used to solve the 
energy equation for temperature. The relative merits of cold and 
hot water injection in oil recovery have been studied here. In 
these problems velocity is not prescribed independently and must 
be obtained from the simultaneous solutions of pressure and 
temperature equations. The following conclusions have been 



arrived at in this study. 


i ) 


ii ) 


iii) 


Raising the formation temperature 
recovery. 

Hot water injection is better than 
injection in terms of oil recovery. 
Loss of heat to the surrounding rocks 


improves 


the cold 


adversely 


the o i 


water 


affect 


the oil production. 
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NOMENCLATURE 


This terminology is applicable everywhere in this work 
except in Chapter 4 and Appendices BCb), E, F and G. 
Nomenclature of these sections is given in Chapter 4, Section 

4.3. 

Symbol s 

x,y Global ordinate and abscissa 

Ax, Ay Grid size along x and y axes respectively 

u,v Fluid velocities along x and y axes respectively 

N,M Total number of nodes along x and y axes respectively 

X,Y Domain dimensions along x and y axes respectively 

t T ime 

At Size of a time step 

Pe Peclet number 

Pe Grid Peclet number 

g 

r) , ? Coordinates parallel and normal to the mean flow 

respectively 
Superscript, 

p Present time step 

~ Vector form 

Subscript 

i,j Nodal indices along x and y axeUltirilsoectively . 
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CHAPTER 1 


INTRODUCTION 


1.1 In-troductlon 

Heat and mass transfer, fluid flow, chemical reaction and 
other related processes occur in engineering equipment, in the 
natural environment, and in living organisms. That these 
processes play a vital role can be observed in a great variety of 
practical situations. Nearly all methods of power production 
involve fluid flow and heat transfer as essential processes. The 
same processes govern the heating and ai rcondi tioning of the 


buildings • 

Major segments 

of the 

chemi cal 

and metallurgical 

industries 

use components 

such as 

furnaces , 

heat exchangers. 

condensers 

, and reactors. 

where thermofluid 

processes are at 


work. Forces that sustain motion of aircraft and rockets arise 
from fluid flow, heat transfer, and chemical reaction. In the 
design of electrical machinery and electronic circuits, heat 
transfer is often a limiting factor. The spread of pollution of 
the natural environment is governed by heat and mass transfer 
principles. Hence there is a need to model and understand the 
processes of heat and mass transfer in the greatest possible 
detai 1 . 

The prediction of performance of a given physical system 
consists of finding the values of the relevant variables such as 
velocity, pressure, temperature and concentration of the relevant 
chemical species in the processes of interest at various instants 
of time. The predictions should also describe the effect of 
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changes in geometry, flow rates, temperature and fluid properties 
on the performance of the system. Heat transfer and fluid flow 
processes can be studied either from laboratory experiments or by 
theoretical modelling. Experimental studies involve tests 
conducted on a full scale or scaled down model of the physical 
system. However, these tests are in most cases expensive and 
time consuming. Further the small scale models do not always 
simulate all the features of the full scale physical situation; 
frequently, important features such as combustion or boiling are 
omitted from the model tests. 

A theoretical. calculation uses a mathematical model 
consisting of a set of differential equations that describe the 
physical system at every point inside the material domain and at 
every point in time. Very often, an analytical solution of these 
equations for a practical problem can not be found, thereby 
limiting its use for an engineer. Fortunately, the development 
of grid based numerical methods such as the finite element method 
and the finite difference method and the spectacular progress in 
recent years in the speed of digital computers hold the promise 
that the mathematical model of almost any practical problem can 
be worked out. This new methodology for attacking the complex 
problems in fluid mechanics and heat transfer has become known as 
Computational Fluid Dynamics (CFD). The main advantage of such a 
theoretical approach are low cost, remarkable speed with which 
the numerical investigation can be performed, complete 
information of all relevant variables and flexibility to simulate 


realistic conditions of our choice. 
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Convection-diffusion equations, arising in heat and mass 
transfer problems, is an important class of problems that is 
being addressed to in CFD. They make appearance in many diverse 
forms such as the IMavier— Stokes equations, Vorticity Transport 
equations. Energy Balance equations and the Chemical Transport 
equation. In general, the transport equations for the velocity 
components, diffusion coefficients and various scalars like 
vorticity, mass fraction and specific enthalpy in a general 
unsteady multidimensional flow are examples of 
convection— diffusion type equations. There are many numerical 
techniques such as upwinding, central difference method and 
hybrid differencing which are available to solve such equations. 
Of all these techniques, upwinding (or upstream differencing) has 
become popular. 

It is well known that the application of the first order 
upwind difference scheme to the convection term of the 
convection— diffusion equation stabilizes computations but 
introduces errors known as numerical diffusion. Owing to this 
error, the solution of the equation shows larger presence of 
diffusion than what is present in the flow. In processes where 
the fluid velocity is high, convection dominates diffusion and 
front like solutions can be produced. In such problems, tracking 
of the fronts is important. Application of upwind technique in 
these problems introduces such a large artificial diffusion 
error that it weakens the effects of convection and leads to the 
smearing of sharp fronts. Fronts can form in certain two phase 
flow problems as well, even when the average velocity involved is 
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small. Resolution of these fronts by upwinding can lead to 
substantial errors unless a highly refined grid is used. 

An example of formation of fronts is seen when the hot water 
injection technique is used to extract oil from an underground 
oil reservoir. The hot water passes on the thermal energy to 
oil, thereby reducing its viscosity and improving its mobility. 
Water does not mix with oil and the former acts as a piston that 
displaces oil trapped in the pores of rocks. This leads to a 
moving interface or front between oil and water. It is 
important that the numerical simulator of the process correctly 
predicts the front movement as the breakdown of front will 
severely affect the theoretical recovery efficiency. The 
nathematical model for the displacement process consists of a 
system of non-linear partial differential equations whose 
solutions exhibit moving fronts. The use of upwind differencing 
tas retarded research in this area due to a host of problems. 
There are the fudging of the fronts, unreliable value of 
/ariables and excessive computation costs of the higher order 
methods that could substitute upwinding. 

There are several other examples where the use of upwind 
technique has not yielded satisfactory results. In a variety of 
•ngineering systems, flows take place in a complex geometry, as 
.n diffusers, nozzles and in corrugated channels, or occur in a 
lanner leading to a complex flow field (as in flow behind a 
:ylinder). Conventional numerical solution of such problems 
ising upwind schemes grossly average out the flow structure even 
it a moderate Reynolds number, say Re 100. This leads to an 



under prediction of pressure drops in such devices. 

Several efforts have been made so far by different 
investigators to come up with a method to reduce upwind errors, 
though only with limited success. This work attempts to find an 
alternative to upwinding itself. A method called the operator 
splitting technique, which completely eliminates the need for 
using the upwind principle, has been used here. Preliminary 
results presented in this thesis show that this new method is 
well suited for both convection and diffusion dominated problems. 

1.2 Llbe Fature Survey 

A review of literature shows that the classical first order 
upwind scheme was first used by Richard Courant. It has been 
explained in detail by Roache C14D. Spalding C13I1 proposed a 
hybrid scheme as an improvement over the upwind scheme. Patankar 
C13!l has recommended a power-law scheme for convective terms in 
flow and heat transfer problems. Hybrid and power law schemes 
derive their inspiration from the exponential solution of one 
dimensional steady flow. This scheme sets the diffusion effect 
equal to zero at much larger grid Peclet number 2) and 
therefore is reduced to the upwind scheme. Leonard C11, C23 has 
strongly criticized the upwind scheme and its derivatives as 
being stable but very inaccurate because of artificial diffusion. 
He has commented that the conclusions of the steady one 
dimensional problem used by upwind, hybrid and power law schemes 
cannot be generalized to more complex problems involving 
convection and any combination of streamwise diffusion, cross 



stream transport in two dimensions, unsteady boundary conditions 
and steady or unsteady source terms. Leonard shows that grid 
refinement can alleviate the problem of artificial diffusion but 
the necessary degree of refinement is often impractical for 
engineering applications. His methods (QUICK and QUICKEST) C23 
solve the differential equation using a third order upwind 
method. Though his method is capable of handling large grid 
Peclet numbers, it is only conditionally stable. In one of his 
latest papers C11I], Leonard discussed higher order schemes such 
as ULTRT— SHARP which is an alternative for high resolution 
non— oscillatory multidimensional steady state high speed 
convective modelling. Murphy C73 applied C* cubic Hermite 
polynomials embedded in an orthogonal collocation scheme to the 
spatial discretization of the unsteady nonlinear Burger equation 
as a model of the equation of fluid mechanics. A new finite 
difference scheme of order accuracy Tree of artificial 
diffusion for solving the steady state incompressible 
Navier— Stokes equations is presented by Her et al . C8I1. The 
scheme uses a false transient approach with a combination of 
variable time step size and over— relaxation Tor the convection 
diffusion terms. Lai et al. C91 have developed a second order 
upwind differencing method for convection— diffusion type 
equations in porous media. The method utilizes explicit 
monotonized upwind/central differencing and operator splitting. 

The operator splitting technique, which is being presented 
here as an alternative to the upwind method was Tirst proposed by 
Yanenko C153. Various possible types of operator splittings 
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feasible in Gas dynamics equations have been described in this 
work. More recently, Issa C3I1 describes a noniterative method 
for handling the coupling of the implicitly discretized time 
dependent fluid flow equations. The method is based on the use 
of pressure and velocity as dependent variables and is hence 
applicable to both the compressible and incompressible versions 
of the transport equations. The main feature of the technique is 
the splitting of the solution process into a series of steps 
whereby operations on pressures are decoupled from those on 
velocity at each step, with split sets of equations being 
amenable to solution by standard techniques. Cooke C43 has 
mathematically established that a certain operator splitting of 
the two dimensional, conservative form, Mavier-Stokes equations 
is second order accurate not only for linear systems, but also is 
true in the presence of non-linearity. An approximate 
(linearised) Riemann solver is presented for the solution of the 
Euler equations of gas dynamics in three dimensions by Glaister 
C57. The scheme incorporates operator splitting and is applied 
to the problem of Mach 3 flow past a forward facing step for some 
specimen equations of state. Ding and Lit C6I] describe a new 
operator splitting algorithm for the two-dimensional 
convection— dispersion-reaction equation. The governing equation 
is split into three successive initial value problems: a pure 

convection problem, a pure dispersion problem and a pure reaction 
problem. For the pure convection problem, solutions are found by 
the method of characteristics. For the pure dispersion problem, 
time— expl i cit finite element algorithm is employed. Analytical 
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solutions are obtained for the pure reaction problem. Results 
are presented mainly for one dimensional problems. 


1.3 Scope of the Present Work 


In the first part 

of 

this 

wor k , 

simple 

one and two 

dimensional test problems 

have 

been 

solved 

using 

the operator 


splitting technique and its results are compared with those of 
upwind and central difference methods and analytical solution 
wherever possible. 

The test problems considered are as follows: 

i) One dimensional transient convection diffusion problem. 

ii) Two dimensional transient convection diffusion problem in a 

rectangular domain with a unit velocity along x axis. 

A temperature step is given at the inlet and the 

transient development of mixing layer is predicted. 

iii) Two dimensional steady convection diffusion problem with a 

unit velocity at 45° to the Cartesian grid. As in (ii), 
two streams of different temperatures flow along the 

diagonal on its either side. The spatial development 

of mixing layer is predicted here. 

The second part of the thesis concerns the problem of oil 
recovery by the hot water injection method. The mathematical 
modelling of the problem consists of three highly non-linear 
coupled partial differential equations, namely two pressure and 
one energy equation. The pressure equations are solved 
simultaneously using finite differences. Fluid velocities are 
computed from the pressure field, which are later used to compute 



the coefficients of energy equation. The energy equation is 
solved using the operator splitting technique. Results have been 
presented to understand the following: 
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(a) Effect of the formation (matrix) temperature 
isothermal (cold water) case. 

(b) Comparison of the hot water and the cold water 
techniques. 

(c) Effect of heat loss to the surrounding rocks on 
water injection efficiency. 


in the 


injection 


the hot 



CHAPTER 2 

DESCRIPTION OF THE OPERATOR SPLITTING ALGORITf^ 
AND THE ONE DIMENSIONAL PROBLEM 


2.1.1 E>escripHon of the Operator Splitting Algorithm 

The operator splitting algorithm (to be called OS) is a 
special case of splitting methods or fractional step methods C15T 
which reduces the solution of a complicated problem to a 
successive solution of simpler problems. Various kinds of 
splittings are possible. These are: 

(a) Geometrical splitting which reduces a multidimensional 
problem to a temporal sequence of problems of smaller 
dimensions or, in particular to transient one dimensional 
problems. In the latter case we shall speak of splitting 
in terms of spatial direction. A.D.I. (alternating 
direction implicit) scheme for the solution of the 
transient two dimensional heat equation is an example of 
this kind. 

(b) Physical splitting in which the original physical process 
is represented as a temporal sequence of processes 
possessing simpler physical structures. This kind of 
splitting can be treated as splitting in terms of physical 
processes . 

(c) Analytical splitting which allows various analytical 
problems to be sol^ved in fractional steps. For example, 
the restoration of divergence in predictor-corrector 
schemes or treating the algebraic terms of the original 
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equations as a separate operator in curvilinear 
coordinates . 

The O.S. algorithm belongs to the splitting of type (C). 
This algorithm properly identifies the mixed mathematical 
character of the governing differential equation and solves each 
homogeneous component of the equation as accurately as possible. 
Hence the OS algorithm solves the governing equation true to its 
character . 

A typical convection-diffusion equation arising in heat 
transfer problems is of the form. 


dT 


A# p 

‘:u. T = ^ T 


( 2 . 1 ) 


In Cartesian coordinates, ^ 


(i- . t->. 


The 
-> 00, 


'TH ’ dy 

mathematical character of the equation is elliptic if t — 

Ar Ar 

parabolic if u is small and hyperbolic if u is large in 
magnitude. In all other cases, the equation is said to be mixed 
in character. These special cases are summarized below! 


A# A# 2 

Steady state (t > oo) (u . '7)T = « V T (elliptic) 

2 

Conduction limit (|u| >0) T^ = oi V T (parabolic) 

j\. A# Ar 

Convection limit (|u| > oo) T^ + u . ^T = 0 (hyperbolic) (2.4) 


( 2 . 2 ) 


(2.3) 


The OS algorithm treats this problem as follows. At any 
time level t, solve Equation (2.4) over a time step At. This is 
considered as a predictor step. The corrector step (Equation 
(2.3)) is solved over the same time interval At. At high values 
of velocity, the correction required is quite small and at low 
velocities, the predictor step acting alone is inadequate. 



However, the two steps combined will furnish the complete 
solution of the Equation (2.1) irrespective of the magnitude of 


u . 


2.1.2 Remarks 

1. The above approach is completely free of upwinding and 
hence devoid of false diffusion errors that usually occur at high 
Peclet numbers. 

E. THe only error arising in the O.S. algorithm is due to 
time discretization and it is of order At. If either the 
predictor or the corrector equation is solved numerically on a 
grid of size Ax, then additional truncation errors Ax*^, n ^ 1) 
are possible. Hence the overall error is, 

e ~ 0 (Ax", At) , n 2: 1 

where n depends on the scheme used to solve the diffusive 
corrector step. Numerical techniques are well established for 
solving the diffusion equation (2.3). This equation does not 
have an analytical solution when solved in a complex geometry. 
Hence for generality a second order finite difference scheme has 
been used to solve for the corrector step. The discretization 
error is then given as, 

e 0 (At, (Ax)^) 

3. The main source of error in solving heat transfer 
problems occur while approximating the hyperbolic terms (u.^T). 
In the O.S. algorithm, these terms are isolated from the complete 
equation and solved as accurately as possible. In this thesis, 
this part of the problem is solved analytically on a local grid. 



4. On non-dimensional ization, the heat transfer equation 

becomes 

T. + (u. ^)T = 7 T (E.5) 

t Pe 

where Pe = UL/« is called as Peclet number <See Appendix B(a)). 
Pe is a measure of the relative importance of convective 
transport with respect to diffusive (or conductive) transport. 
Most of the existing algorithms fail or become inaccurate when Pe 
is raised. The OS algorithm proposed here is uniformly valid 
over any range of Pe. Its application to engineering problems 
where Pe can be small or large forms the focus of this thesis. 

As the Peclet number is raised, a discontinuity in the 
initial conditions will propagate through the flow domain 
unchanged. This is called^ front. Predicting the occurrence of a 
front and its location is a challanging numerical problem. 

Fronts can also form in non-linear problems at low or moderate 

Peclet numbers. The oil recovery problem studied in this thesis 
is an example of this (See Appendix B (b)). 

5. The OS algorithm is only conditionally stable since it 
involves explicit time marching and is limited by the Coura n t 
stability principle. The stability criterion is stated in words 
as , 

"the numerical movement of the front should not exceed one 

grid spacing in one time increment since the maximum front speed 

is the fluid velocity". 

At 

Hence, we require -j— 


i 1, for stability. 
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(In dimensional form, the above condition is, — r— i 1) 

Ax 

See Appendix C for the derivation of stability criteria 
using von-Neumann analysis. 

2. 2. 1 The One PimensiQnal Problem 

The application of the OS algorithm for a one dimensional 

problem is given below. 

The governing differential equation is 

T., + T = -J- T (2.6) 

t X Pe XX 

subjected to initial condition: t = O, T = 0 and the boundary 

conditions : 

X = 0, T = 1 

X > oo, T > O. 

This represents a physical problem in which a hot fluid with 
unit velocity is introduced suddenly (at t = 0) at the inlet of a 
flow domain €See Figure2.1>. The fluid in the physical domain 

(0 < X < oD) is initially at a constant cold temperature. If the 
Peclet number is large, a thermal front whose thickness is small 
will move through the flow region with unit velocity. As 

Pe > 0, the front will be diffuse and the energy transport is 

due to the molecular action of thermal diffusivity. 

2. 2. 2 Analytical Solution 

The closed— form solution of the above equation is. 


00 



T(x,t) 


•Trr 


-Cx?16T)^>Pe> ^ 
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The above integral is well behaved and is evaluated using 
Simpson's rule for numerical integration. Values of the upper 
limit between 8 and 10 reduce the integrand to an order 10 ^ for 
any Peclet number and x studied here. 


2.2.3 Central Difference Formulation 

The governing differential equation is discretized by a 
fully implicit method in time as. 


tP+1 



At 


T. .) 
1+1 1—1 


p+1 


SAx 


1 xJi + l ^i-1 

' ;i;:7 


the 


The above equation determines the value of temperature 
i^*^ node. fSee Fig. S.1>. It can be rewritten as, 


at 


ft xr’ * B T 
1-1 


p*'' . c rPi:; = D T 
1 1+1 1 


where 


A = 


1 

2Ax 


Pe(Ax) 


B 


1 

At 


Pe{Ax) 


11 1 

C = - 2 and D = — 

2Ax Pe(Ax) At 

For a grid system shown in Fig. 2.1, the matrix structure 
arising froam the central difference/implicit formulation given 
above is as follows: 
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The node N is located far enough for the condition = 0 to 
hold. The above matrix is inverted using a tridiagonal algorithm 
( TDMA ) . 


2 . 2 . 4 Upwind Formulation 

To generate diagonal dominance, 3.t is essential to use the 
upwind scheme for the convection term, T^. This is as follows. 


u T 

X 


u T. - u T. ^ 
e 1 w 1-1 

Ax 

u - u T. 

e 1+1 w 1 

Ax 


, u > 0 


, u < 0 


Subscripts *e* and *w’ denote the walls of the control 
volume which are to the east and to the west of the node i 
respectively. In the present problem u = 1 (>0) and therefore 
the finite difference formulation becomes. 


-P+1 


At 


jP+1 _ yP+1 

i ‘ i — 1 

Ax 


Pe 


, ^i+1 ^ ^i+1 ^^i ^p+1 

( ) 

(Ax> 


We write this equation as, 
{ « T._, B T, + C T, 


= D T? 


for each node i and 
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invert 


the matrix 
A = 

B = 

C = 


using TDMA. Here, 

— - 1 
Ax Pe(Ax) 

11 2 
+ + j_ 

At Ax Pe(Ax) 

1 

^ and D = 

Pe(Ax) 


1 

At 


Since |Bj > jAj + |Cj, the matrix generated by upwinding has 
strict diagonal dominance when Dirichlet conditions are 

prescribed at x = O and x > oo. The inequality becomes 

stronger as At is reduced. 

2.2.S Operator Splitting Algor ithm 

The Equation (2.6) is split as. 

Predictor: + T =0 (2.7) 

t X 

and 

Corrector: T^ = T^^^ (2.8) 

with each equation applied over the same interval At. The 
boundary conditions are relevant only over the second corrector 
step. The convective step Equation (2.7) has an analytical 
solution, 

T(t-x) = constant. This can be written as, 

T (t + At - X) = T (t - (X - At)) 

Choosing At = Ax, this becomes, 

((t + At) - X) = T^ (t - (x - Ax)) 
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Hen 

ce. 

the 

temperature at a 

point at 

the 

new t 

ime 

step 

(p+1 

)is 

obta 

ined 

by transferring 

the valu 

e prevailing 

at 

the 

cur r 

ent 

time 

(p) 

at the rearward 

node. If 

At 

and Ax 

do 

not 

mate 

h, then 

the 

value of T in the 

interval 

(i-1. 

i > at a 

dist 

ance 

At f 

rom 

node 

i i 

s transferred to 

i • 






The 

cor 

rect 

or step is solved 

by finite 

differences 

using 

an 


implicit time marching scheme. Hence, 

- T? 1 T._^. + T. . - ST. p+l 

1 1 _ ^ 1 + 1 1-1 

At Pe (Ax)^ 

or 

+ B T. + C T. = D T.’ 

1+1 1 1-1 1 

where 

IS 11 

A = - . , B = — + C D = — 

Pe(Ax) At Pe<Ax) Pe(Ax) At 

Once again TDMA is used to invert the matrix arising from 
the above equation. Successive applications of the two steps 
completes the calculation for a single time increment At. 

2 . 2.6 Results 

Figures 2.S, 2.3, 2.4 and 2.5 show the plot of temperature 
distribution as a function of distance for various time levels 
and Peclet numbers. Typical values of Ax and At used in these 
calculations are 0.1 and 0.1 respectively. All four approaches 
given above are compared in these figures. The agreement among 
the numerical predictions and the analytical result is very good 
at low Peclet numbers (Pe £ 10). At Pe = 100, only the operator 



splitting algorithm matches the analytical result closely- The 
upwind method fails due to false diffusion inherent in its 
formulation- The central difference method fails due to matrix 
ill conditioning and loss of diagonal dominance at high Peclet 
numbers. The operator splitting algorithm is seen to be 
uniformly accurate at low as well as high Peclet numbers. 




One dimensional problem 









Figure 2.5 : One dimensional problem : Pe =100 



CHAPTER 3 


STUDY OF THE TWO DIMENSIONAL PROBLEM 

3. i • 1 Eteficr Ip-Lion 

We consider a two dimensional convective heat transfer 
problem in which the flow is given as uniform. The governing 
differential equation for this problem is, 

’’’t "’’x = PF~ ^^xx 

with the boundary conditions 

0<x<X, y = o ,T = 0; 

0 ^ X < X , y = Y ^ j ^ ^ . 

X = O , 0 :< y < Y/2 , T = 0 j 

X = 0 , Y/2 < y < Y , T = 1; 

x = X .OiyiY ,|I=0 

and the initial condition 

0<xs:x, 0:£y<Y, T - 0. 

See Figure 3.1 for the description of this problem and the finite 

difference grid. Two parallel streams of equal unit velocity 

(u = 1,v s 0) but unequal temperatures (T=0 and T=1) are taken to 

come in contact suddenly at the inlet of the rectangular domain. 

Simultaneously, the upper wall (y = Y) is maintained at a 

temperature T=1. The fluid in the physical domain is initially 

at a constant cold temperature T=0, possessing the aforementioned 

unit x-direction velocity. A mixing layer will form at the 

interface of the two streams in which the temperature gradually 

changes from the higher value to the lower one. The thickness of 

this layer grows in the downstream direction. At high Peclet 
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numbers, the mixing layer is thin and the temperature 

distribution is almost a discontinuity. If Pe > 0, the 

mixing layer spreads quickly to touch the boundaries. 

In this problem, the flow domain taken to be square whose 
edge is 2 units. We have studied temperature profiles in the 
mixing layer at x = 0.32 and x = 0.82 at a time t = 0.5. In 
these calculations a 26x26 grid is used with Ax = Ay = 0.08 and 
At =0.1. The grid size is chosen so as to adequately resolve 
the front in any given problem. 


3.1.2 Upwind Formulation 

The discretization of the convection term is along the 
same lines as that of the one dimensional problem which was 
discussed in section 2.2.4, Chapter 2. This is repeated he^e. 


u T 

X 


u 


e 


T. 

1 


u 

w 


T 


i-1 


Ax 


u 

e 



u 

w 


T. 

1 


u > 0 

u < 0 


Since u = 1, therefore the finite difference formulation 
becomes. 


jP+'l _ jP 


^p+1 _ .|.p+1 


T. .. .+T, ^ ,-2T, 


i , j i ,j ^ irj _ f i + l y j i~1 , j i y j 

At Ax Pe t (Ax)^ 


T- + T. . . 

i,J+1 i,J-1 


- 2 T. 
1 


(Ay) 


p+1 

-} 


(3.2) 
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We write this equation as, 

■IaT. ..+BT. .+CT. , .+D +ET \ 

\ 1-1, J i,j 1 + 1, j i,j-1 i,j + 'IJ 


P+1 


= tP 


1 , J 


Here , 
A = 


Ax 


At 


Pe(Ax) 


2 ' 


B 


for each node i 


1 H. , + 2^ 

Ax Pe(Ax>^ Pe(Ay) 


At 


Pe(Ax) 


D = E = 


At 


Pe(Ay)‘ 


We have used the Gauss— Seidel method to invert the system of 
equations generated by Equation (3.2). Notice that the 
Scarborough criterion of diagonal dominance is unconditionally 
satisf ied, i . e. , 

|Bl> |A| + |C| + |D| I- |E| . 

Hence the Gauss— Seidel iterations converge unconditionally 
in this problem. For any node (i,j), the iterations are 
generated by the formula. 


T , 
1 , J 


1 r J 


- (A T. 


i-lx j 


+ C T. . . + D T. . . + E T. 

1 + 1tJ i,J-1 i,J + 1 


B 


Here, the superscript denotes the previous iteration 

value and the new value. 

For the right hand side boundary, where the Neumann boundary 
condition is imposed, we have 


This is 


T 


i+1, j 


incorporated 


= T. . at i = N. 
1 , J 

in Eq. 3.2 to yield. 
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T. . 
1 r J 


- (A T. ^ + D T „ + E T. J 

1 > J 1-1 ,J I , J-1 I > J 1 

B + C 


3.1 .3 C^>«!rat,or SpliHing Formulation CADi;> 

The Equation (3.1) is split as. 

Predictor: T* + T =0 (I) and 

t X 

Corrector: ^ (T^^ + T^^) (ID , 

each equation applied over the same interval At. The boundary 
conditions are relevant only over the second corrector step, as 
in the one dimensional problem. 

Solutions of the predictor and the corrector is calculated 
as follows: 

Step 1: The convective operator (I) has an analytical solution, 

T(t-x) = constant. This can be written as, 

T^'^^((t + At) - X) = T*^ (t - (X - At)) 

As in the one dimensional problem, the temperature at a 
point at the new time step is obtained by transferring the value 
prevailing at the current time at the rearward node if At = Ax. 
If At ^ Ax, interpolation is required. 


Step 2: The corrector step consists of the two dimensional heat 
equation. It is solved using the efficient Alternating direction 
implicit (ADI) technique. The ADI method is second order 

r 2 2 21 

accurate with a truncation error of Oj (At) , (Ax) , (Ay) j. 

In ADI, we first compute row by row intermediate 

temperatures while sweeping along the ’j’ axis -fsee Figure 3.1 




Then, by sweeping along the ’i’ axis, final temperatures at the 
end of a timestep are calculated in the column by column 
fashion. TDMA is applicable here for computing row or column 


temperatures. The timestep for each sweep is ~ At/2. 

The discretization for each sweep is given below. 


Sweep along 'J* axis: 


jP+1/2 _ ^p 

i , j i . 


(At) 




ET. . p+1/2 


(A x) 


x)2 j 


f ^i>j+1 ^ ^i,j-1 ^'i,j 1 ] 

I (Av)2 j j 


|a T. . . 

\ 1-1, J 


left i,j 1+1 




p+1/2 


Bright T. .+D T. 

^ i,J i,J+1 


+ E T 




where 


A = C = - 


Pe (Ax) 


®left 


2 (At) 


Pe (Ax ) 


D = E = 


Pe(Ay) 


right 


1 - 


2(At) 


Pe(Ay) 


Sweep along *1* axis: 


.p+1/2 


(At) 


T. . . . + T. 


- 2T. 


r r i+i,j i-i,j i,ji 

U (Ax)^ J 


p+1/2 


f ^i,j+1 ^ ^ItJ-I ^'i,j' ] 1 
t Pe(Av)^ J J 


Pe(Ay) 
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or , 


/at. . . 
\ i,J-1 


®left ^ ^ ^ 


}" ' {®ri 


T. . + D T.^, 
ght i,j 1 + 1, j 


+ E T 


i-1,j } 


p+1/2 


where 


■<At) 


= C = 


adi 


2(At ) 


Pe (Ay ) 


®left = ^ 


Pe(Ay) 


adi 

T~ 


(At) 


D = E = 


adi 

T 


2(At> 


Pe (Ax ) 


' ®right ^ 


adi 

T 


Pe (Ax ) 


In the OS algorithm, successive applications of the two 
steps completes the calculation of temperature for a single time 
increment At . 


3.1 .<4 The OS Algorithm < I<jaus;s Seidell* 

To check the results predicted by ADI, the diffusive 
corrector step has been solved independently by a Gauss-Seidel 
method. 

The corrector step equation in finite difference form is. 


- tP . 

1 r J 1_»2 

At 


1 f T.^^ .+T. . .-2T. . 

I 1+1, J 1-1, J i,J 

— j ^ 

Pe ^ (Ax) 


^i,J+1 ^ ^i,j-1 ^^i,j 1 

(Ay)^ / 


2T, , >.p+1 


or. 


/a T. ^ . + B T. . 
\ 1-1, J i,J 


+- C T._,. . + D T. , . + E T 

1+1, J i,J-1 


i, j + l} 


P+1 


= tP . 

1, J 




where 
A = C 

B 



-At 

? 

Pe (Ay ) 

2(At ) 

Pe(Ay)^ 


Here as in upwind formulation, 

|B| > I«1 + |C| + |Bl |E| 

Matrix inversion proceeds along the same lines as in the 
upwind method (Section 3.1.E). 


3. 1 .5 Results 

We have solved the corrector step of the OS algorithm 
using two different numerical techniques to ascertain that the 
difference between the predictions by the upwind and the OS 
algorithms is not due to the use of different matrix inventers 
such as the Gauss— Seidel and the ADI methods. There is 
practically no difference between the three solutions namely OS 
(Gauss-Seidel ) , OS (ADI) and upwind at Pe = 0.1, 1 and 10 -^^See 

Figures 3.2 and Figure 3.3^ for both cross sections (x = 0.32 and 

X = 0.72). At Pe = 100 and Pe = 1000, we do observe a difference 
between the OS and the upwind predictions. Here, the OS (ADI) 
profile continues to be identical to the OS (Gauss— Seidel ) 
profile. The temperatures predicted by the upwind method are 
higher than those obtained from the OS algorithm. This 
difference can be attributed to longitudinal false diffusion 


arising from upwinding. 
















*0.72, t*0.5 
0 S. (Gauss-Seidel) 
O.S. ( A.D.I.) 
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Figure 3.3 ; Two dimensional problen 


at X 


0.72. 


Pe = iO 



_J i_ 

1.28 1.60 



_Pe^OOO 

■0 B 8 0 5 8 ir 


1.28 


= 0.1 


i 


: ir-e 


1 / 1 







3.2.1 Grid-Orientation Effect 


In the one and two dimensional problems studied here, we 
observe the effect of false diffusion in high Pe flows which are 
aligned with the Cartesian grid. In addition to the severe 
truncation error in computations involving field variable 
curvature in the streamwise direction, there is also the problem 
of streamline to grid skewness in two and three dimensions C1,E3. 

This is called the Grid-orientation effect. Although this 
problem can be alleviated by skew differencing procedures C12I1, 
numerical experience indicates that the remaining truncation 
error is of the same order as the skewness error that has been 
eliminated. In addition, there is a possibility of convective 
instability, which makes it inadvisable to use the technique in 
modelling the convection of transport variables. Thus, the 
relative cost— benefit ratio of skew differencing is rather 
unfavourable . 

An approximate expression for the false diffusion 
coefficient of the upwind scheme for a two dimensional problem 
has been given by deVahl— Davis and Mallinson (197E); it is 


false 


p U Ax Ay sin 20 
_ — 

4 (Ay sin 6 + Ax cos 


C131 


where U is the resultant velocity, and 0 is the angle between O 
and 90°) made by the velocity vector with the x direction. See 
Appendix A for the role of ^^^jse causing truncation error in 
the upwind scheme. 
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o 

Note that _ is a maximum when & - 43 

false 

of skewness, upwinding is expected to be least 
demonstrate below that the OS algorithm continues 
under the circumstances. 

Consider the two dimensional heat transfer 
the velocity field is not necessarily parallel to 
governing equation is. 


u T + V T = 
X y 


Pe 


(T + T ) 
XX yy 


Me specify 


At this level 
accurate. We 
to be accurate 

problem where 
the grid. The 

(3.3) 




so that the flow is at 45° to the grid. 


The boundary conditions are as shown in Figure 3.4. The 
initial condition 0<x<X, 0<y<Y, T {x,y,0) = 0 corresponds 
to an initially cold fluid. 


This 

problem is identical to 

the 

two 

dimensional 

problem 

described 

in Section 3.1.1 except 

for 

the 

fact 

that 

the two 

parallel 

streams of equal unit 

velocity 

and 

of 

unequal 


temperature (T = 0 and T = 1)» instead of being aligned with the 
grid, are inclined to it at 45°. Initial temperature of the fluid 
in the physical domain is T = 0 (except at the boundaries where 
Dirichlet boundary conditions are prescribed). Here again, a 
thermal mixing layer is formed. For definiteness, we study the 
problem at steady state. At steady state, the mixing layer 
thickness is inversely related to the magnitude of the Peclet 


number. 
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Steady state is reached at a time level t = 2.8 units. The 
following data is used in the numerical computation: 

At = 0.04, Ax = 0.04, Ay = 0.04 

Grid size - 51x51 
Region size — ExE 

3 . 2.2 Altai yllcal Solution 

For an infinite region, the steady state temperature 
profile can be obtained at high Peclet numbers by solving the 
reduced equation. 


where r) and f are 
flow. 

The boundary 
7) = 0 , 
r) = 0 , 
as 

as 

The analytical solution of this problem is, 

T(r),? ) = (1 + erf §) 

where 

? Pe 

* 2 J Y) 


(3.4) 


V Pe 

coordinates parallel and normal to the mean 


conditions are, 

? > 0 , T = 1 

? < 0 , T = 0 

? > 00 , T > 1 

? > - <», T > 0 


This is valid for Pe >> 1. 
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3.2.3 Upwind Foriraila'Lion 

Discretization of Equation 3.3, using the upwind method for 
the convection terms, results in the following nodal equation: 


.p+1 




1 > J 1 r J 

At 


1 


r T. . - T. . . 

[ 

V Ax 


T. . - T. . p+1 


1 » J 


Ay 


} 


1 

Pe 




T. . ..+T, ,^p+1 

(Ay)^ 


This simplifies to. 


•f A T . . . 

\ 1-1, J 


+ B T. . + C T.^. . + D T. . . + E 

IrJ 1 + 1, J 1,J-1 


i, j + l} 


P+1 P 

= T. . 
1 , J 


for each node i, where 


At 


A = 


At 


At 


At 


T2<Ax) Pe(Ax)' 


, B = 1 + 


|2(Ax) 

2(At) 


At -At 

2 ' ^ ~ 

Pe(Ax) 


Pe(Ax) 

At 


4H(Ay) 

2(At) 

Pe(Ay>^ 

At 


42 (Ay) Pe(Ay)' 


, E = - 


Pe<Ay) 


The Gauss-Seidel method is used to solve for . at all 

1 jr J 

the internal nodes- 


3w 3. 4 Op erator Splitting Formula-tion CADI!> 

The governing equation (3.3) is split as, 

1 1 

(I) T^ + -r=- T + T =0 (predictor) and 

t 42x 42y 
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(II) 


1 

Pe 


(T 


XX 


+ T ) 


yy 


(corrector) . 


The analytical solution of step I is obtained by using a new 
1 

coordinate 4 > = ^ ^ (x + y) . 

This result in. 


T^ + T. = O whose solution is T(t -<?!>) = constant, 
t 4 > 

This can be written as 

<t + At - <(6) = T^ (t - (^ - At)). 

For Ax = Ay and At = A<^, this means that the temperature at 
a point at the new timestep is obtained by transferring the value 
prevailing at the current time at the diagonally opposite 
upstream node, i.e.. 


tP-"'' 

1 72 


tP 


1 » j- 


1 


The corrector step II is solved as in Section 3.1.3 using 
A.D.I. 

3.2.5 Operator ^littlngf Formulat ion COauss-SeidelD 

Here step II is solved by the Gauss-Seidel technique 
instead of ADI. See Section 3.1.4 for details. 


3. 2. 6 Results 

Comparison of profiles obtained by the analytical solution, 
upwinding, OS (ADI) and OS(Gauss— Seidel ) is shown in Figure 3.5. 
The analytical solution cannot be used for comparison at low 
Peclet numbers such as Pe = 10 and Pe = 100 due to 

a) the inapplicability of the boundary-layer approximation 
and b) the finite size of the region considered in numerical 
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simulation. At Pe = 1000, the mixing layer is thin and 
temperature reaches values 1 and O well within the flow domain. 

At Pe = 10, both upwinding and operator splitting give 
identical results. The temperature profile becomes sharper at Pe 
= 100 and this is well reproduced by the OS solution. For Pe = 
1000, we find that OS results come close to the analytical 
result. Note that at Pe = 100 and Pe = 1000, the upwind results 
are almost identical. We conclude that for the grid used, false 
diffusion dominates over physical diffusion in the upwind 
solutions at Pe > 100. 

We also notice that OS (ADI) and OS (Gauss-Seidel ) profiles 
are almost identical in the three cases. Hence we are certain 
that the differences in the O.S. and the upwind results are due 
to false diffusion rather than differences in the Gauss— Seidel 
and ADI methods. 

For the fully evolved mixing layer, it can be shown that the 
mixing layer thickness varies as 

<5 4 Pe = constant. 

In Table 1, the 6 -fPe values as predicted by the two 
schemes are given for three Peclet numbers. We see that upwind 
scheme values come nowhere close to being a constant. OS scheme 
shows much better response vis— a— vis the analytical solution. 

In Table 2, S 4 Pe values predicted by OS method for 
different grid sizes is listed. We observe that despite a change 
in the magnitude of the constant, <54^ remains nearly constant. 
This means that an increase in Pe reduces <5 in the right 
proportion and produces a suitably sharp temperature profile. 



43 


Tabl« 1 X Comparison of <5 -fPa Valuas 


Pe 

Analytical 

Upwind 

OS 

10 

13.01 

6.86 

6.66 

100 

13.08 

19.32 

7.81 

1000 

12.88 

30.173 

8.95 

Tabl«» 2 

X Effaet of Grid 

»-coar Sorting on 

Vai u«ps 

Pe 

N(=M) = 51 

N(=M) = 41 

N(=M) = 31 

Ax(=Ay) = 0.04 

Ax(=Ay) = 0. 

05 Ax(=^y) = 0.067 

10 

6.66 

6.64 

6.614 

100 

7.81 

7.954 

8.014 

1000 

8.95 

7.826 

10.14 



CHAPTER 4- 


NUMERICAL MODELING OF HOT-WATER INJECTION METHOD 

IN OIL RECOVERY 


4. 1 Int/r-oducLlon 

Before we start with the application of operator splitting 
algorithm in the numerical simulation of the hot-water 
injection method, a brief description of various oil recovery 
processes is first presented C163. A common misconception about 
hydrocarbon recovery is that oil or gas is found in large pcols 
in underground caverns and must be pumped out of the reservoir in 
a manner similar to pumping a liquid out of a storage tank. This 
is not the case. In general, the hydrocarbon is trapped in the 
microscopic pores of a rock, like sandstone, and will flow 
through the rock only under the influence of extremely large 
pressure differentials. A large percentage of pores are 
connected and the fluids can flow through these linked pore 
paths- However, the paths are very small, irregular and twisted. 

In many instances, the pressure found in reservoirs are 
extremely high. Often these pressures are high enough to force 
the resident fluids to flow through the porous medium and out of 
the wells without any pumping effort at the wells. As the 
reservoir pressure is depleted through fluid extraction, the 
hydrocarbon production decreases and finally stops. This type of 
recovery, termed primary recovery, extracts 15— 30/i of 


hydrocarbons in the reservor. 
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Subsequently, a fluid such as water may be injected into 
some wells in the reservoir while petroleum is produced through 
others. This serves the dual purpose of maintaining high 
reservoir pressures and flow rates and also of flooding the 
porous medium to physically displace some of the oil and push it 
towards the production wells. This type of pressure maintenance 
and waterflooding is usually termed secondary recovery. 

Unfortunately, waterflooding is still not extremely 

effective and significant amounts, upto 50X, of hydrocarbon often 
remain in the reservior. Due to strong surface tension effects, 
large amounts of oil are trapped in small pores with narrow 
throats and cannot be washed out with routine waterflooding 

techniques . 

In order to recover more of the residual hydrocarbons, 

several enhanced recovery techniques involving complex chemical 
and thermal effects have been developed. These techniques form 
part of a variety of methods termed tertiary recovery and 

have been given the name enhanced oil recovery CEOR) techniques. 
Since oil is trapped in the pores of the reservoir because of 
surface tension, several EOR techniques involve the lowering of 
the surface tension to allow the oil to flow from the small 
pores. Some of these techniques are given below. 

1) Surfactants flooding: This process involves the use of a 

surface active chemical which gets into the pores and causes 
reduction in the interfacial tension between oil and water. 
This way capillary forces are reduced and it leads to 

improved oil production. 

2) Miscible flooding: It involve injection of gases such as COp 
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and other chemical solavents which mix with the resident 
hydrocarbon. Together they form one fluid phase of lower 
surface tension. Once the mixing or "miscibility** is 
attained, the fluids will flow together in one phase, 
eliminating distinction between phases, and enhanced 
hydrocarbon recovery is possible. 

3) in situ combustion: When the hydrocarbon reserves are too 

deep to be mined, techniques for in situ transformation of 
the hydrocarbons into states which can be pumped to the 
surface via production wells are required. In this 

technique, the solid hydrocarbon is ignited in the ground and 
oxygen or air is pumped into the injection wells to maintain 
combustion. The? hot hydrocarbon gases produced in the 
combustion process are pumped to the surface via the 

production wells for use as low-grade hydrocarbon. 

4) Thermal recovery methods: This techniques for lowering the 
viscosity of heavy oils involves addition of very high 
temperatures to effectively "melt" the viscous hydrocarbons, 
allowing them to flow more readily through the porous medium. 
A primary technique for thermal flooding is to inject either 
hot water or steam into the reservoi"* via certain injection 
wells. As the heated fluids come into contact with the 
cooler oil, the temperature is raised. Consequently the oil 
viscosity is reduced and it can be pushed easily through the 
porous media. Getting the high temperature fluid down into 
the reservoir is a difficult engineering problem. The hot 
injecting fluid may lose heat while moving through the 
injection well and later to the host rock bounding the oil 



rich region. Some energy is also consumed in the chemical 
reactions of hydrocarbons. Apart from heat loss, there is 
also the problem of gravity override. This occurs when the 
hot fluid rises on to the top of the oil saturated sand layer 
and bypasses it. If the flow rate is high, the interface 
between the resident petroleum and the invading water can 
become unstable leading to the formation of long fingers. 
These fingers grow in length towards the production wells, 
^yP^ssing much of the oil. These factors cause a drastic 
reduction in oil production. 

It is this last type of EOR method based on thermal recovery 
that we are concerned with in this chapter. We numerically study 
EOR using hot water in a small two dimensional rectangular test 
cell. The operator splitting technique is used to solve the 
energy equation. 

<1. 2 Literature Review 

Hot water injection technique is a tested technique used in 
petroleum industry. Many attempts have been made in the past to 
model the fluid and heat flow involved in this problem. 
Researchers have taken recourse to analytical methods to predict 
the movement of saturation and temperature fronts in the oil 
fields. Lauwerier*s E21I] attempt to formulate the equation of 
heat balance and determine its analytical solution for a simple 
real life situation is one of the earliest efforts in this 
direction. Spillette LZOl has refined Lauwerier’s approach and 
extended it to steam injection as well. There are a variety of 
theoretical methods listed in References C16I3 and C17D, but they 
are essentially extensions of Lauwerier’s and Spillelte’s works. 
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Despite the elegance and apparent simplicity of analytical 
solutionsT they are plagued with problems. Analytical solutions 
are possible only for very simple physical situations. When real 
life difficulties such as heterogenet ies in the soil, multiphase 
multicomponent fluid flow, distillation effect of hydrocarbons 
are taken into account, no analytical solutions are possible. At 
best, the analytical approach can yield initial estimates of the 
parameters in the physical processes. Numerical modelling of the 
reservoir process is an answer to the limitations of the 
analytical approach. The computational approach allows us to 
include as much physics as possible of the processes involved. A 
large number of coupled nonlinear partial differential equations 
can be solved simultaneously by numerical methods. Heterogenet ies 
and geometrical complexities present in any reservoir can be 
taken into account. Moreover it produces detailed information 
about the process which is useful to the oil engineers. As a 
result of these advantages, numerical modeling has rapidly gained 
importance . 

Reference C19I1 gives an excellent introduction to the 
general topic of reservoir simulation. Description of various 
oil recovery methods is given by Ewing C16I3. It is also a good 
source of references of important research papers in this area. 
Boberg C17I1 describes thermal methods of oil recovery in detail 
and presents an extensive compilation of case studies and field 
data. He also presents two state— of— art numerical simulators 
CH23 and CS33 for steam injection process. An important paper by 
Stevenson et al C24I] describes some of the latest computational 


methods in reservoir simulation. 
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4.3 Nomenclature 

x,y Global grid coordinates m 

Ax Ay Grid size along x and y axes respectively m 

X,Y Overall domain dimensions along 

X and y axes respectively m 

X .,Y Problem domain dimensions along 
d d 

X and y axes respectively m 
i,j Node indices along x and y axes respectively 

^ •A. 

i, j Unit vectors along x and y axes respectively 

t Time s 
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V 

oy 

V 

wx 


V 

wy 


T 


T 

g 

U 

V 


D^rcy velocity of oils y component m/s 
Darcy velocity of waters x component m/s 
Darcy velocity of waters y~compoanent m/s 
Temperature 
Formation temperature 
Generation temperature at source 
Composite velocitys x direction m/s 

Composite velocitys y direction m/s 
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4. 4 Def 1 nl t i ons 

Here we explain some of the terms used in this chapter. 

(a) Porositys It is the ratio of total pore volume to the total 
volume of porous medium (sand). 

(b) Percentage pore volume (ppv)i It is an index of efficiency 
of oil recovery, ppv is the ratio of total volume of oil 
produced to the total pore volume. It is calculated as 


ppv = 


Initial volume of oil— Present volume of oil 


X 100 


Pore volume 


nd M 


1 , ) 

j_.| ' initial 'present 

HH R 

2 2 e 

i=1 j=1 


X 100 


nd M 

y y (S <i,j)| . _ . , - S (i,j)| . ) 

Z. Z O I initial o I present 

i=1 j=1 

(nd . M) 


X 100 


where S is oil saturation, 
o 

(c) Oil and water velocities: v and v are the volume 

Q W 

averaged velocities called Darcy velocities and are related 
to the local-fluid velocities as follows: 

V = Se V, V =S£-V. 

o O O W WW 


(d) Saturation: It is the fraction of the space available to 
flow (i.e. pore volume) occupied by a phase. 

j rr ^ 
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(e) Biot number (Bi): It is an index of heat loss to the host 

rock bounding the oil rich sand ^ See Figure 4.1 Higher 

the value of Bi, higher the heat loss- 

(f) Oil displacement efficiency: is the efficiency with which 
oil is flushed out from the porous medium by water. It is 
directly related to ppv. Higher the ppv, higher the output 
of oil and consequently higher the oil displacement 
efficiency. 


(g) 

Overhang: 

It is the impermeable rock lying on 

the 

top of 


the porous 

medium -j^See Figure 4.1^ and is 

also 

called 


"overburden 

i". 



(h) 

Underhang: 

impermeable rock lying underneath 

the 

porous 


medium and 

is also called "under burden". 



( i ) 

Formation: 

refers to the whole structure consisting 

of host 


rock and oil rich sand. 

4.5 Physical Description 

- The physical domain of the two dimensional problem is shown 
in the Figure 4.1. Essentially, we would be studying the water 
piston effect in a rectangular domain packed with sand and 
guarded by impermeable but thermally conductive rocks above and 
below the porous region (called overhang and underhang 

rerpectively) . Initially the porous medium (i-e. sand) is 
impregnated with oil to a given saturation. It has some initial 
temperature (called as formation temperature) and is at rest. 

We assume in this study that the pore volume not containing oil 
has water. At equilibrium, the oil and water pressures in sand 
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a.r6 uni'Quely clet^6rmin©d “froro t^h 0 ir s^tlLursit^ion . Oil r 0 cov©ry is 
initialBcl by pushing w^ter at a given higher pressure and a given 
temperature We observe the movement of saturation and 

temperature profiles in the porous region. These profiles would 
have the shape of fronts because the Peclet number 

in this problem is quite high |see Appendix B<b) 

The following two variants of the main physical phenomena are 
addressed to in this study. 

(i) Isothermal Problems Here the water that is being pushed in 
is of the same temperature as that of the formation (T^). 
Therefore, we need not solve the energy equation. 

(ii) Non— isothermal Problems Heat is transported into the flow 

domain through hot water that is pushed in. Here, the 

temperature distribution is important and the 

energy equation must be solved along with the pressure 
equations . 

4. 6 Scope of bhls Chapter 

Results have been presented in this chapter for the 
following cases. 

(i) Effects of formation temperature recover/ 

( ppv) . 

(ii) Comparison of oil recovery for the isothermal and 
non-isothermal cases for a given value of T^. 

(iii) Evolution of saturation, temperature and velocity profiles 
for the isothermal and non-isothermal cases. 




( i V > 


Effect of heat loss, to the overhang and the underhang, on 
oil recovery. 
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- 4 . 7 Governing Differential Equat ions 

Since the phenomena of two phase flow through the porous 
media is too complicated to be described on a microscopic scale 
equations applicable on a macroscopic scale have been used in the 
present study -^See Appendix D for a discussion on the 

macroscopic approach^. These equations are derived in Appendix 
E.The pressure equations for oil and water phases are. 


Oil: - 


S ^ 
o o St 




dS ^ p Sp 


S ? — - 

o ^ o ^t 


w 


w 


‘Ip, 


St St 


ow 


k k p 

V . ( — — —) 


^ P, 


(4.1 ) 


Water: 


-S ft ^ 
w ‘ w St 


^P, 


w w St 


w 


dS Sp Sp 


w 


w 


dp. 


St St 


ow 


w 


k k p 
rw w 

7. ( V p 


(4.2) 


Equations 4. land 4.2 are strongly nonlinear coupled 
differential equations. Their form is that of a parabol 
conduction equation with source terms 


S ft -yj— 

W W «t 


dS 


w 


dp. 


.1^ 

^t 


Sp 


w. 


'ow 


partial 
ic heat 


and 
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Here 


Tempe ra-lure is determined -from the energy equation. 


SJ 

sx 


+ u 




£I 

Sy 


< Li_) 

‘"t 


(4.3) 


U = 


V p c 
ox o o 


V p c 
wx w w 


|(1-S ) P c + S P cl 
^ w oo w w wj 


+ (1-£)(pc), 


V = 


V pc + V pc 

oy o Q wy w w 


-|(1-S ) p c + S pc 1 
w o o w w w J 


+ (1-€)(pc) 


R 


The following equations are also valid. 


p = (S ) = p 

” f.i ” 


ow 


w 


w 


(4.4) 


S + S =1 

o w 


(4.S) 


Velocities are determined using Darcy’s law as. 


w 


k k 


ro 


V = ~ 


V = ~ 


k k 


rw 


w 


(7 p^) 


(V p ) 


(4.6) 


The equations of state are. 


P_ = P_ 

o o 


F ref "{ 


•I ‘Oo - fo.ref’ - ^ef) 


} 


(4.7> 
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’w - '=w,rsf { ' * <P„ - P„,ref> - - T,,,>} 


^o,w 


k = k (S ) 

ro,w rOrW w 


Pc " Pc ‘V 

OW OW 


See Table 4 for the dependence of k , p on T, 

ro,w 

OW 

S , S respectively, 
w w ^ 


4. B Iriitial and Boundary Conditions 

The initial and boundary conditions, along with much of 


thedata ^ Table 3^- is adapted from Boberg C171, p.24. These 

values have been used by Boberg for a laboratory test where the 
dimensions of the test cell were of the same order as the domain 
dimensions of this problem. 


At the left boundary of the domain 4 See Figure 4 






water 


at a pressure of 190 psi and at a temperature of 100 C is 
applied. Saturation of water is 0.86 here. |^It is not unity 

because of some residual oil which can not be flushed out^. It 


is the maximum water saturation possible anywhere in the domain- 
At the right hand boundary of the porous region, water pressure, 
temperature and saturation are maintained at 190 psi, 50 C and 
0-2 respectively. 

Initially, the porous region is maintained at a uniform 
water pressure and temperature of 190 psi and 50 C respectively. 
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Qil saturation is O A nil 

IS obtained by adding water 

pressure to the capillary Dres<;iirck tKa t 

F xciry pressure, the latter being a unique 

function of water saturation. 

The initial conditions of p^, p^, and T were smoothened 
out to avoid the computational instabilities resulting from a 

step like initial condition. These conditions are given below. 

At time t = 0 , 

= (260 - 175 X) psi; Q < x < 0.1 X 


- 190 psi 

= 0.86 - 1.65x 

= O.Z 

T(x,y) = (100-125x)°C 

= 50 °C 


0.1X < X < X 
0 < X < 0.1 X 
0.1 X < X < X 

0 i X s: 0.1 X 


0.1 X < X < X 


At time t > 0, the impermeability conditions at the bounding 
surfaces of the host rock give. 


'^oy '^wy y ~ 0,Y 

Using Darcy’s law -fEquation 4.6, Section 4. 


the above 


impermeability condition is stated in terms of pressure as. 


at y = 0, Y. 


Further, at x = O, 


S^(0,y) 


= 260 psi (= 1.793 MPa) 


= 0.86 


T (0,y) = 


100°C 


Po(0,y) 


p^(0,y) + p^ (S^(0,y)) 
ow 
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At X = X, 


P„<X, y) 

= 190 psi 

(= 1.31 MPa) 

y> 

= 0.2 


T(X,y) 

II 

Ul 

o 

o 

n 


p^(X,y) 

= Pw(X,y> 

+ p^ (S^ (X,y)) 


ow 


Boundary conditions related to heat loss to overhang and 
underhang is given by: 


ft 

Sy 


+ Bi. (T - T^) 


O at y = O, Y. 


4r. 9 Discretization of P ressure Equations 

The oil pressure equation I^Section 4.7i Equation 4.1^ can be 


written ass 

dS dp dS 

-S/9^4-(S? <^Ji - 

o o o o dp St dp 

. k k p 

= _L vp 


^p, 

dt 


w 


The discretization of equation is fully implicit in time, 
i.e., all spatial derivative are evaluated at the new time level. 


k k p 
ro o 


^ • ^o'^w' • ^o'Po,''' 

e.p„(T) 


Let e 
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1 

+ 

ifj 


u- 


p -p P+1 

O 0 

i,j + 1 


p -p 

o. . o. 


P+1 


Ay 




A y 


Here &, & r & t ^ are evaluated as harmonic averages of 
e w n s 

their values at nodes that lie on their side (Figure 4.1). 

The left side of the pressure equation is discretized as. 


d S 


(S ? 
o o 


w 


dp. 


) . . 
1 , J 


ow 



Po""'' - '■o 1 

IrJ IrJ 

-f- 

r ds "1 

w 

■ 

At 


dp^ 




*'ovir 


^p+1 

w . . 

IrJ 


w . 


It J 


At 


^ S ft 
o o 


1 T J 
Jp .JP-I 

i,j i»j 


At 


The full equation in discretized form is written. aS: 


I A p + B p + 

o o .... o *^0 . . 

i + 1 fJ i»J 


C p + D p Pn 

O ^W. . □ i ° 

i,J 1-1, J 


i, j + 1 
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where 


+ F p 
o 


i f j-1 


p+1 


= G 


(4.8) 


A 


r ® "1 

® 1 

li 

o 

a 

Cl 

A 

dS 1 

w 

p 1 

(Ax)2 J 

dPc 

^ OW 

. . ' At ’ 

1 T J 


dS p 

gP V _ ( ) 

-*0 dp^ ^i,j 


1 , J 


B. 


ow 


At 


It J 


e + e Q + e '\ 

e w_ ^ _n s 


(Ax) 


(Ay)‘ 


D 


-f 


e 1 ^ 

r 

^ 1 

w E - - 

I 


'"o. . J ° I 


(A y)2 J 


1 r J 


1,J 



Backward differencing in time of the temperature term is 
Justified because the main effect of T on pressure distribution 
is through fluid properties. These do not change as rapidly 
with time as the primary variables S. 
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The water pressure equation (section 4.7, Equat 
also discretized as: 






ion 4.2> is 


\ w Pin + D p + E p 

L '•'i-Hio « '-i^j Oi,j w wPw.^ 

P+1 


J+1 


+ F 


w 


i » j-1 


= G 


w 


(4.9) 


where 

A 

w 





4 c = J 

r ‘^^w 

1 




. (Ax)^ 

i _ *1 

J " 1 

L ‘“’c 

^ OW 

* At 

— 


D 


w 


- { r^)' ■■ ■ 

'■ w. . ^ 

Ir J 


i-J '• ■^>-4 


At 


w . 

IwJ 


{ y + ^ 

^ e ^ w 


n s 


>P 


(Ay)^ 


W 


(■ 


i r j 




(Ay)' 

i r j 


M- 


S ? 

w . w 

IwJ 


dS 

r i!L 1 

I dP J . . 

^ C ^ IrJ 
OW 


. dS . 

[-FFT-li.J V. 

OW 


At 


w . 


IrJ 
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w . 


1 tJ 


w 


JP Jp 

i>j 

At 


where 


k k p 
rw 



k^ (S ) .p (p ,T) 

rw w 

(T) 


The discretized version of the boundary conditions at y=0, Y 
(See Section 4.Q) are. 


p (i,1) -p (i,2) = O 

o o 

P„ '‘-1> - P„<‘-2> = O 

(i,M— 1> — p^(i,M> = 0 

p - p (i,M) = O 

w w 

Dirichlet boundary conditions at x = 0, X are easily applied by 
replacing Equalions(4.8,46) by the respective function values of 
the pressure. 


-4. 10 Solution of the Pressure E^uat-lons 

When the equations 4'6 4-3 of Sec. 4.9 are written for 
each node and such pairs of equations are drawn up for all the 
nodes in the domain -fsee Fig. 4.lT a set of linear equations is 


obtained. The matrix structure is obtained: 
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The band width of the matrix is and its size is 
(2 NM) X (2 NM) . This matrix is inverted using a sparse matrix 
solver C29D that requires values of only those el ements that lie 
within the band. 

^.11 Solution Cft blue Enorgy Equabion 


The energy 

equation 

is solved 

using the 

operator 

splitting 

algorithm. 






The energy 

equation 


= (-ii> 9^1 



at 


. v|I = 

ay 

is split 

into the 

following steps. 






Predictor ! 


Sx 




Correctors 

at 

t ^ > ^'^T 

"'t 

» 



as described in 

Chapters 

£ and 3. 





4r. 1 1 • 1 Solublon of bhe Predicbor using Sbreamlinos 

We choose a curvilinear coordinate system — Y) ) locally 
at each node such that is aligned with the net fluid composite 
velocity at that node, i-e. ? is the streamline passing through 
the node |see Figure 4. 19, (A) | f and T) are orthogonal to each 

other. Then, it can be shown that 

U T + V T = U' T„ where U* = 

X y s 



See Appendix F(a) 



(xi,yi) 
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Hence, the predictor step becomes, 

• . . 

+ U =0 where U = U 

The solution of this equation is, 

V 

T (U t - ) = constant 

Since the grid surrounding each node is 
streamlines ? to be st raight , i . e . , 


<x,y) . 


small , 


we 


assume the 


= x* cos<^ + y’ sin(^ |^See Figure 4.19, (D)| 
-1 V 


U 


and - tan 

Here (x*,y’) are the local coordinates at (i,j) and (x,y> are the 
global coordinates for the grid. 

The solution of the predictor step requires that the 
temperature at *0* at p^^ timestep be transferred to the point 


P* at the (p+1)^^ timestep -{See Figure 4,. 19 




(D)j., 


1 • e . 


T^"^^ (U’ (t+At) -?) = T^ (U’t - -A^)) 


where A^ = At. 

- See Appendix G for the details of interpolation required to 
determine T^. 

The predictor step is implemented along the following 
lines, 

(a) From U and V, calculate ^ and U’ = -I 

(b) From the value of find out the quadrant in which Q lies 
tc> Xq and y^ are calculated using the following relations 

A^f = U* At 

A Xq = A? . COS^, Ay^ = A^r .siiW^ 

Xq = Xp - AXq , ~ 
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(d) Tq is computed by interpolation and transferred to Tp 

(e) Steps (a) — (d) are repeated for every internal node in the 
flow region. 

N. B. At the boundaries, y = 0 and y = L, V s 0 

4. 1 1 . 2 Solution of the Corrector Step using A. D. I. 

The corrector equation is, 

} 

This equation is solved using ADI as described in Section 

3.13 (Chapter 3). The only difference is that unlike Pe in 
Kh 

Section 3.13, ( ) is nota constant and therefore it needs to be 

"^T 

calculated at each node. 

4.12 Algorithm 

The system of equations governing the distribution of oil 
and water pressures and temperatures are non-linear and mutually 
coupled. They must be solved simultaneously and by iteration. 

The algorithm, in short, is as follows! 

1. We start with the initial p^, p^, T and S^ values at t = 0. 

2. The coefficients (A — G and A - G ) of pressure 

o o w w 

equations are computed using the values of p^, p^^, T and at 

the latest timestep. 

3. A system of pressure equations, as shown in Section 4.10* is 
solved to obtain the new values of and p^. 

S , p and p are updated using the new pressure values, 
w "^o w ^ 


= (- 


2 

•)7 T 


from Section 4.11 


4 . 
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5. 

Darcy 

velocities v , 
ox 

V , V 

oy wx 

» V are 

wy 

calculated using 

the 

latest 

values of p , p 

I S , p 

W *^0 

and p . £See 
w 

Appendix G> . 

6. 

Using 

S r p T p,, , V , 

wow o 


coef f icients 

of the energy 


equation are computed and the equation is solved using the OS 
algorithm yielding a new temperature T. 

7. and p^ are updated once again using the new T. 

8. Steps S to 7 are repeated until convergence of p , p and T 

■^o w 

is achieved. 

9. Fresh computation is initiated for the next time level 
starting from step 2. 

■4.13 Results arK i Conclusion 

The data for oil and water properties and those of the 
porous formation used in this problem is given in Tables 3 and 4. 
This data is adapted from the data used by Boberg C17!l, for a 
laboratory test which reproduced steam and hot water injection 
processes on a small scale. 

The domain which is used for calculating ppv and for 
plotting the scalars has a length (= 1.6 m) which is smaller 
than the overall domain length X (= 4 m) -^See Fig. 4.l|-. In 

this manner, the influence of the right hand boundary conditions 
is delayed. On a smaller domain, in the isothermal case, there 
is a build up of oscillations in the saturation profile. For a 
total domain size of 4 m, oscillations do not appear for a period 
of 1 hour of simulation and the flushing effect of water on oil 
is seen over a distance of 1m. 

The computer program developed in the course of the present 
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study showed the following performance. On Convex Minisuper 

computer, the isothermal version of the program took two minutes 

of CPU time for a 41x5 grid and a time step At = 0.01 hour. 

The non-isothermal version took three CPU minutes. Convergence 

P+1 P P+1 

criteria used was (Q - Q )/Q x 100 < 0.01 where Q is d . d 

^o ^w 

and T. The convergence was achieved, on an average, after 
five iterations. In order to minimise the computational 
instabilities, smoothened initial conditions were used (See 
Section 4.8)- Further, timestep was gradually increased from 0.1 
(At) to At, where At is the characteristic timestep. 

The correctness of the computer code and the matrix inverted 
has been tested by solving the linear transient heat conduct:ion 
equation . 

Now, we would be presenting the results in a format given in 
Section 4.6. 

(i) Isothermal problem: Figure 4.2 shows the rise in saturation 
level of water in the flow domain along its centreline with tinte. 
There is no front and the profile is diffusion dominated. Next, 
we consider the effect of raising the formation temperature. As 
we raise from 30°C to 50°Cy we find an increase in the water 
saturation along the centreline of the flow domain (See Figi.re 
4.3). Me also observe in Figure 4.4, that oil recovery ifriproves 
with the raising of the formation temperature. This is expected 
because of the drop in oil viscosity with an increase in 
temperature. Evolutions of centreline velocity profiles of water 
and oil are shown in Figures 4.5 and 4.6 respectively. 

(ii) Oil displacement efficiencies of the isothermal (cold water 
injection) and non-isothermal (hot water injection) cases are 




Figure 4.2 : Isothermal case : Progress of saturation front 
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gure 4.6 : Isothermal case i Change in oil velocity pro 
with time* 





SS9Jd 


Figure 4,7 : Non- iso thermal case : pressure profiles 
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compared in Figure 4.8. The percentage pore volume (ppv) of the 
non—isothermal case is below that of the isothermal case for 
times less than 1 hour. Beyond one hour, the performance of the 
non-isothermal process is better than that of the isothermal 
process. A possible explanation for this phenomena is given in 
part (iii) below. 

(iii) We analyze next the hot water injection technique in 
detail. Figures 4.9 and 4.10 show the development of centreline 
saturation and temperature profiles. Both of them show a 
distinct front like structure and on comparing them, we observe 
that saturation front moves with the same speed as the 
temperature front. The saturation profile reveals an unusual 
feature namely a dip in immediately beyond the front. This 
dip is responsible for the initial low level of the ppv of the 
non-isothermal case as compared to the isothermal case fSee 
Figure 4.8>. 

Figures 4.11 and 4. IS give the development of water and oil 


velocity profiles respectively for the non— isothermal case. 


We 


N 


find that the magnitude of velocities is higher than that of the 
isothermal case. The profiles individually show the presence of 
a peak which travels with the same speed as the two 
aforementioned fronts. 

(iv) We now study the effect of heal loss to overhang and 
underhang on oil recovery. As mentioned earlier, (Section4.5) 
this is effected through the manipulation of the parameter Bi . 
Higher the value of Bi, greater the heat loss. Since this is 
convection dominated transport phenomena (See Appendix B(b)), 



isothermal 




001 



Figure 4,9 : Non isothermal case : Progress of saturation fron 
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Figure 4.11 i Non isothermal case : Changes 

profile with time . 
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for the time frame considered (" 1 hour), the effect of heat loss 
does not penetrate deep enough to influence the centreline 
characteristics significantly for the domain width considered (Y 
= 0.8m). At smaller widths <Y = 0.04 m. See table 3) 

the influence of heat loss is seen to be significant. 

Figure 4.13 shows the effect of increase of Bi on 
centreline saturation profile. As the heat loss increases (i.e. 
Bi increases), the saturation front profile becomes increasingly 
smooth. Study of the centreline temperature (Figure 4.14) 
reveals the cause of the decrease. As Bi is raised, a smaller 
amount of thermal energy is available causing the weakening of 
the centreline temperature front. This leads to a higher 
averageA^iscosity and so a lower fluid mobility. Hence the volume 
of water being pushed through the sand is reduced, and 
consequently, water saturation profile is lowered. 

In Figure 4.15, the typical widthwise temperature profiles 
for various Bi values are shown. When Bi = 0.01 and 1, the 

profile is a straight line which implies that the overhang and 
the underhang are perfect insulators. But as Bi is raised to 100 
and then to 10, 000, the profile increasingly becomes parabolic in 
nature. The slope is a maximum at the top and the bottom of the 
domain because of the heat loss to the rocks. 

The effect of heat loss on ppv and the oil displacement 
efficiency is shown in Figure 4.16. We observe that for t < 1.5 
hour, oil recovery (ppv) for Bi = 0.01, 1, 100 10,000 are 

nearly equal. At the end of 2 hours, a clear trend can be 
discerned. We find that oil recovery is the least for Bi 





hour 
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Figure 4.16 ; affect of Bi on oil displacement efficiency 




10,000 and a maximum for Bi = 0.01 and 1. Hence we conclude that 

in the long run, heat loss from the sand bank during hot water 

injection adversely affects the oil production. 

We now discuss the effect of the size of the grid on 

resolution of fronts of various scalars. In Figure 4.17, the 

consequence of progressive grid coarsening (41x5, 31x5 and 21x5) 

on the centreline profile is studied. The difference in the 

profiles is striking for the cases 31x5 and 21x5, but the 

differences become smaller as we move on to a 41x5 grid. We 

observe that it becomes increasingly difficult to keep track of 

the saturation front faithfully with grid coarsening. Not only 

does the front lose its true shape, its speed also becomes 

erratic. Same conclusions can be derived from the study of 

centreline temperature profiles (Figure 4.18). Results 

presented in this chapter use a 41x5 grid. 

Finally, the role of the field data in the form of K (S ) 

ro w 

and K (S ) and p (S ) (Table 4) must be mentioned. Here, 

rw w '^c„ w 

ow 

.. given for a finite set of S values. A linear 

ro,w ^ w 

interpolation is used to compute intermediate values of the 
relative permeability and capillary pressure. Recent literature 
(Reference C25I]) shows that the prediction of a numerical 
simulator for flow through a porous media is sensitive to the 
values of K , K and p . Any deviation from the correct value 

these inputs can lead to unrealistic predictions. It is 
suggested here that closed form analytical expression be used to 
furnish the relative permeability and capillary pressure values 
as in Reference C273. These are expected to improve the quality 












08 
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grid size on tem; 

: non-isothermal ca 
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gf the numerical predictions. 

Before closing this chapter, we would like to mention the 
^vantages of harmonic averaging over the arithmetic one of B and 
1..0C (See Section 4,9). The reason behind the choice of 

Y S/OlI ^ 

averaae is that viscosity has very sensitive exponential 

harfnonic. 

dependence on temperature and the relative permeability changes 
sharply with saturation. This can cause sharp variations in B 
and values in a small region. The choice proved to be a 
correct one as the wiggles in saturation profiles in the 
isothermal case were reduced in magnitude or not present at all. 
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Table 3 : Data Used in the Problem : Parameter Values 

Most of the values in part (B) are taken from pE24, Boberg 

C171. 

(A) Domain length, X = 4.0 m 

Domain width, Y = 0.8 m 

DoMin width. Y = 0.04 m (To study Bi effect) 

Domain grid size. Ax = 0.1 m 

Domain grid size. Ay = 0.2 m 

Domain grid site. Ay = 0.01 n, (To study Bi effect) 

Domain of analysis ('problem domain' : Figure 4.1). 

X, = 1.6 m 

d 

time step. At = 36 s 

(B) Permeability, k = 132 Darcies 

Porosity, ^ = 0.375 

S'* . n-T JF = 5x10"^ 1/psi C=0. 03447 1/Pa3 

Compressibility* Oil, ^ 

compressibility: water.(;„= 3.1x10'^ 1/psi i:=0.0Z137^1/Pa3 

Expansivity: Oil. » 0.00041 1/°F C=Z.28x10^^ 1/°C1 

+ = o 00041 1/^F E«28x10 1/ 

Expansivity! water, 0- 

P = 0.5 Btu/lbm-°F C=H092.0 

Specific heats oilr 

W-s/kg- 

tar c = ^ Btu/lbm-°F C= 4184.0 

Specific heat: water, 

W-s/kg 

Heat capacity/volume of the porous medium. ^ ^ 

(pc)^ . 36 Btu/ft^-°F C = 241EW6.S H-s/m - C3 

tha Dorous medium. 

Thermal conductivity of t P ^ 

,s* /ft C= 0.1661 W/m- Cl 

K =40 Btu/ft-day r i- « 

•^h 
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Injection temperature, T^= 100°C 

Formation temperature, = 50°C 
Injection pressure (water) = 260 ps 
In-situ pressure (water) = 190 ps 


i 

i 


C= 17.926x10^ PaD 
C= 13.1x10^ Pa3 
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Table 4 ’ Data Used in the Problem ; 

K CS3, p CS3 Tables 
ro w c w 


(i) ' 


T(°F) 

^^(cp) 

50 

5000 

150 

1Z6 

S50 

20 

350 

7.1 

450 

4.0 


(i) P^<T) 


T(°F) 

IJ ( cp) 

'^w 

60 

1.130 

75 

0.935 

80 

0.875 

100 

0.685 

120 

0.560 

150 

0.430 

200 

0.308 

250 

0.230 

300 

0. 182 

400 

0.145 

500 

0.120 


CS 

W 
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0. 1 

0 

1.0 

4.111 

0.2 

0.0016 

0.875 

0.095 

0.3 

0.0081 

0.735 

0.072 

0.4 

0.0259 

0.590 

0.061 

0.5 

0.0672 

0.420 

0.051 

0.6 

0.1000 

0.210 

0.041 

0.7 

0.1400 

0.070 

0.031 

0.8 

0.2000 

0.016 

0.021 

0.86 

0.2500 

0 

0.011 


For intermediate values, linear interpolation was used to 

calculate oil and water viscosities and oil and water 

relative permeabilities ^'^ro *^rw^‘ While computing 

saturation from capillary pressure p , for p^ ^ 0.95 psi, 

ow ow 

exponential interpolation was used. For p^ < 0.95, linear 

ow 

interpolation was used. 



CHAPTER 5 

CONCLUSIONS AND RECOMNENDATIONS 


The performance of operator splitting algorithm was tested 
against those of the upwind and central difference methods by 
solving simple one and two dimensional benchmark problems. 
Results were compared against the analytical solution of 
problems which were available in some cases. It was found that 
for low Peclet number flows, the performance of OS was identical 
to the performance of the other two schemes and it matched the 
analytical solution (whereever it was available). For high Pe 
cases only the results by the OS came close to the analytical 
solution. Upwinding showed egregious +'alse diffusion errors and 
central differencing was also not satisfactory. Hence we have 
concluded that the new scheme is relati;vely free from the basic 
errors afflicting the conventional lower order schemes namely 
false diffusion and grid orientation errors. OS was also 

successful in tracking fronts in high Pe flows. 

The new technique was successfully used in modelling hot 
water injection technique for enhanced oil recovery and the 
effects of various parameters on oil recovery were demonstrated. 
They were as follows: 

improvement in oil production by raising the formation 
temperature T^ . 

- superiority of hot water injection method over the cold water 
approach . 

adverse effect on oil production of heat loss to the 


surrounding rocks. 


There is one feature of the new algorithm which needs 
improvement. We have solved the hyperbolic part of the 
convection-diffusion equation analytically using the formula 

T(x, t + At) = T(x — At, t). If At = Ax, we obtain by 

ix) ^ 

p 

equating it to interpolat ion between 

and is used. Interpolation on a coarse mesh causes 

errors which are similar to false diffusion {See Figure 5.1>. 
As shown in the figure smaller the value of At/Ax , higher the 
error. A numerical scheme could have been used to solve the 
hyperbolic part, but the advantages of accuracy of the 
analytical solution would then have been lost- Figure 5.1 shows 
that interpolation <?rrors diminish on a refined mesh. 

Scope tar Future Work 

Before declaring the operator— spl itting method as an 
efficient alternative to the upwind scheme, it needs to be 
tested on various other flow configurations. Some of these are 
convection— diffusion equations with source term and radioactive 
decay and to fluid flow problem. In the context of oil 
recovery, the OS algorithm should be compared with those which 
are already in use in the industry. 
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APPENDIX A 

FALSE DIFFUSION 

Consider a one dimensional convection diffusion equation. 






■) where r = diffusion 


coefficient and its solution on grid is shown as: 




On discretizing using Central difference scheme for the 


convection term we get 


_ 




^E- 


+ rf ^ 1 

I A 2 J J 


Where u, r and Ax are constant, ^ is the average value of 
the cctntrol volume. 

Now we use the upwind scheme for the convective term. 


** r *C ^ ( *E *U ^C] 

5t = [ — g J 

Ax '■ Ax ^ 


where 


<^i - ifu>0 

L w 

= <^g. if u < 0 


Equation (ii) can be recast into the same form as Equation 
(i> by writing 


“ — 
2Ax Ax 


( i i i ) 
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On comparing (iii) 
diffusion coefficient due 
given by the formula 

r 

num 


with ( i ) » 
to upstream 

and is 


the additional numerical 
differencing is seen to be 

positive irrespective of 


the sign of u. 

r insignificant in comparison 

ThuSy wishes to mak xium 

r the condition required 
.ith the physical diffusion coefficient r, 

I lA /T This is clearly much more 

n /p = Pe = I u|Ax/I^ KKZm 

'^nufc' 9 ‘ P» < 2 for the 

.trin9ent than the practical stability condition Pe^ - 

Both conditions are infact highly 

central differencing. Both co 

. .tiral grid requirements for problems 

unrealistic in terms of practical g 

- 1 9„tprest in which the global 

of engineering or geophysical interest, 

A nn a tvoical macroscopic leng 
Peclet or Reynolds number based on a yp 

be extremely large. 


scale can 


APPENDIX B 


(a) Non-dimensionalization of the Heat Transfer equation 


The energy equation is of the form. 


^t “ '' ST 


Z Z 
Sx Sy 


Let L = characteristic length and U 
The characteristic time is L/U. 

Let X* = x/L, t*= tU/L, 


= characteristic velocity. 

y* = y/L and 

T - T . )/(T -T . ) 
min max min 


be the nondimensionalized variables. 

Then the energy equation in dimensionless variables is. 


* * 
U ^T U * ^T 

L ^t* L “ dlT* 


U * ^T* 
L '' Sy* 


(V*2 T*) 


Hence 


SJ 

SX 




1 

Pe 





(i) 


where Pe = — 

o( 

The superscript can be dropped in further calculations 

once it is understood that u,v,T, x and y are dimensionless • 

{ b) Non dimensional izat ion of energy equation of Chapter 4 
The energy equation used in Chapter 4 has the form. 
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where U and V are ’composite' velocities constructed from the 

individual water and oil velocities. Let U, L, T be the 

characteristic composite velocity, length and time. T = ^ . 

U 


In dimensionless variables, we get 



u 

* 

JLw* ^ 
* 

L 


L Sy 



7 T 


Hence 


dt 


+ U 


£I 


+ 




1 



* 




where 

Pe 


L U 



The length scale L cannot be properly identified in this problem. 
However, simulations in this study are carried out over a region 
size of 4m and we choose L=1m for the moment. 

For a typical case 

= 2,642, 548.8 , U = 1-65x10 ^m/s 

T 3 o« 

m ~ C 

and so Pe = 2424.8. 


Hence 

we conclude 

that 

Peclet number in 

the 

hot water 

injection 

method is of 

the 

order of 2500 and 

the 

flow 

i s 

convection 

dominated. 

As seen 

in the solutions of 

the 

one 

and 


two dimensional test problems, this high Pe flow will 
sharp fronts. 


generate 
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APPENDIX C 

von-NEUMANN STABILITY ANALYSIS OF THE PREDICTOR STEP IN 

OS ALGORITHM 


The question of weak stability of the predictor would be 
answered using von Neumann or Fourier analysis. See for 

details. 

The predictor step is, tP^^ = tP . - ) (tP .-tP ^ .) 

1 ,J l,j Ax 1,J 1~1,J 


Let the error be e(x,t) = T b <t) e 

L m 


ik X 
m 


Let us study the growth characteristics of a single error 


term e (x, t ) b ( t) e 
m m 


ik 

ro 1 


Now b^ is chosen as e^^ where k = ^^^is the spatial 

m m L ^ 


wave number. 


at 

On substituting (ii) in (i) and cancelling the e^ e 


term, we get 


a ♦ * ♦ A -ik Ax 

aAt , uAt. f m 

e <= G) = ^“"Xx ^ I ^ 




a(t+At) 


where G = 


is the amplification factor. 


Hence G=(1-i.») + v e where v = and /9 = k Ax 

Ax m 


|G| - j ( 1— V ) + V cos/? — ivsin/3| 


J {1-^(1-cos/9)>^+ ii>sin/9>^ 


than 


After a little manipulation, this reduces to: 


lOS 


jQj = jl + (1 - COS/?) .2u . (Vi - 1 > 

For stability, the amplification factor |Glshould be less 
1. Hence (1 - cos/? > .2^ . <v-1 ) < 0. 

Since w > O, we require v < ^. 
the condition for stability is 


O < 


(^) < 

Ax 


1 


Hence 
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APPENDIX D 

THE NEED FOR A MACROSCOPIC APPROACH [183 

In principle, "the continuum equations that describe various 
transport pahenomena are known and may be written at the 
microscopic level within the porous region. Here we focus our 
attention on what happens at a point within a considered phase 
such as oil or water present in the porous soil. We may even 
know the conditions that prevail on the surface that bounds the 
phase. However, at this level, the equations cannot be solved, 
since the geometry of the surface that bounds the phase is either 
not observable or too complex to be described. The same is also 
true for point values of variables within the phase. As a 
consequence, the description and solution of a transport problem 
at the microscopic level is impractical and perhaps impossible. 

In many branches of science, it is convenient to ignore the 
particulate nature of matter and to adopt the hypothesis that 
matter is a hypothetical substance that is continuous throughout 
the spatial domain it occupies and can be described in that 
domain by a set of variables which are continuous and 
differentiable functions of spatial coordinates and time. This 
continuum model of matter serves as the fundamental postulate of 
continuum mechanics and thermodynamics. 

The continuum approach, so useful in treating problems 
related to a single phase, can be extended to a multiphase system 
such as the porous medium, where the various phases are 


separated from each other by abrupt interfaces. Accordingly, the 
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real system, consisting of two, or more phases, that together 
occupy disjoint subdomains within a porous medium can be 
replaced fay a model in which each of the phases is assumed to 
behave as a continuum that fills up the entire domain. This 
domain is termed as the macroscopic space. For each point within 
this macroscopic space, average values of phase variables are 
taken over representative elementary volumes (REV), centred at 
the point regardless of whether, in the real domain this point 
falls within that phase or outside it. The averaged values are 
referred to as macroscopic values of the considered variables. 
By traversing the entire porous medium domain with a moving REV, 
thus assigning averaged values to every point, we obtain fields 
of macroscopic variables which are differentiable functions of 
the space coordinates. 

Apart from the advantage of circumventing the need to find 
exact configuration of the interphase boundaries, the other 
advantages of the macroscopic approach are as follows. 

- It describes processes occurring in porous media in terms 
of differentiable quantities, thus enabling a solution by 
employing methods of mathematical analysis. 

- The above mentioned quantities are measurable. 

These advantages are at the expense of the loss of 
detailed information concerning microscopic configuration of 
interphase boundaries and the actual variation of quantities 
within each phase. But the macroscopic effects of these factors 
3re still retained in the form of coefficients, whose structure 
and relationship to the statistical properties of the void space 


(or phase) configurations can be analyzed and determined. In the 
case of porous media, the numerical values of these coefficients 
must be determined experimentally, in the laboratory, or in the 
field. 

Altogether , in view of the advantages of the continuum (or 
the macroscopic) approach, it shall be employed in Chapter 4 to 
describe transport phenomena in porous media. 
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APPENDIX E 

derivation of the differential equations governing oil-water 

FLOW IN A POROUS MEDIUM 

We derive here the equations governing oil and water 

pressures and temperature in two phase flow through a porous 

medium 

Assumpt 1 OB:ts I 

(i) Two fluid phases flow simultaneously. 

(ii) The fluids (oil and water) are immiscible, i.e. there is 
no mass transfer between the fluids. 

(iii) The energy balance assumes instantaneous thermal 
equilibrium between the fluid and rock matrix. This is 
justified since the phase velocities are small. 

(iv) Gravity term is neglected while calculating Darcy 
velocities. It is reasonable because the densities of 
oil and water are close and the time frames are short. 
As a result, the gravity override effect will not show 
up in this study. 

(v) The surrounding rocks and the matrix are incompressible 

and have constant transport properties. 

(vi) Relative permeabilities of phases is a function of 

saturation and viscosities are dependent on temperature 
only. 

(vii) There is no net mass source in the flow field. 

(viii) Specific heats of the fluids are independent of 


temperature and fluid pressure. 
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(ix) Viscous dissipation term in energy equation is neglected. 


Derivations 

Following are the basic equations which are used to derive 
the final pressure and temperature equations for flow and heat 
transfer in a porous medium. The assumptions given above are 
used in these derivations. The nomenclature for this Appendix is 
the same as that of Chapter 4. 


Darcy* s Laws 


k k 'v 

ro 


(Vp ) 


(i ) 


w 


k k 


rw 


w 


(V p ) 


(ii) 


Mass balances 


Oil 


o o 
St 


+ V . V P = q 
o o 


(iii) 


^(£S p > 

Water : + V . v p q, , 

St w w w 


(iv) 


Since the rock is incompressible/. £ = constant. Further 
and q , the mass sources are taken to be zero in this work. 


Energy Equations 


Oil and water: 
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7 .K^7T 


-w.lv p H + v p Hl=4r rc(S p H -fS o H 
1^0 0 O W W W J «>t L O O O WWW 


+ (1-£)(pc)„T +qcT +q c T 
R ^0 o g ^w w g 


Since H = c T and q = q . s 0, we get 
® o,w o,w o w 


7K^ - V 

n 


.fv p c T-**v pc t]= -y= Tecs p c t+s p. c t) 
o o w w J L ° ^ ° WWW 


+ <1-E) (A>c)p T j 


Equations of state: 


J Sp 

= _JL -Jl-i 

w ^w ^Pw It 


P 


1 dp 

^ |t ' ^o 


1 0 I 

o I 0 


p (S ) = P„ - P.. 

ow 


<vii ) 


( vii i ) 


We encorporate the equations of state into the mass balance 
equation for oil. Using (iii) and (vi), we get 


dp dS„ 

£^0 ^ ^ ^ ■ '^0^0 = ° 


r po ^ 

-j^o P^ dt" dt 


. dS 


— 4- ^ . V p =0 
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Using Darcy’s law (i> in this equation, we get 


6=0 '=’0 

r, - 

^ 1 

Hh 

6 

D 

li 

0 

CO 

< 

k k p 
ro 0 

^ Po 

L 0 ^t 

^0 ax J 


*^0 ax 

^0 








( ix) 

Similarly for water. 






r 

ap 




as 

k k p , 

rw w ^ ^ 

^ w ^w |_ 

^ w ax 

^w ax J 

+ 

£ • 

'^w dt 

^w 

*'w 








(X) 

Now we would be 

deriving 

the 

two final 

pressure 

equations 

which have 

to solved 

simultaneously. 



Rewriting (ix) and (x) as 







a p as 

*^o . 0 


1 L 

k k p 

ro'^o „ 

p1 

“S ft -yr 

0 0 

+ s ? 

0 0 

^t ^t 


e^o 1 
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ap 

^ + 

as 

w 
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^Pw} 

-S ft jrr 

W W 

+ s ? 

w w 

^t 

ax 

£/=>« 1 

^w 

and noting that p 

(S ) = 

w 

Po “ 

■ Pw ' 




ow 






d(1 - S ,) 

0 __ w _ 

d(1-S„) 



1 

T 



ax 

‘‘Pc 







'"ow 








“ ^w 



►Po ^Pw 1 





dp^ 


^t ^t 





ow 






ds 

w 

ds r ^ 

w 

p ^ 

Pw 





^t 

dp j 
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we get 


Oili 


-S /9 ^ 


^Pc 

^0 ^ 


dS 


w 




“OWM 








1 


k k p 
ro o 


P P 


^ p. 


(A) 


Similarly, 


Waters 

— S /> -SY 4- 

w w St 

Sp 
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dS dp 

+ —ii (-12 

dp^ 'dt 

dp 

_ '*') = 

dt 






ow 




1 

= . . 


r k k p 
rw w 

1 ^ n 


(B> 


^W 



J Pw 




(A) and (B) are the final form of pressure equations. Now, 
we will transform the energy equation from the conservative form 
(v) to a nonconservat^ive form using mass balance equations (iii) 
and (iv). This is necessary to use the OS algoriothm for solving 
the energy equation. We have. 


Oil 


r d(sp ) 

£ L_f_ 

^ St 


V . V p 

o o 


= O 


. C T 
o 


f ^(S p ) 
Water : J g w w 




^ . V p = 0 f . C T 

w w I w 


On adding these two equations, we get 


^t 


('I 


< 1 - S ) c (O + c 
w o o 


S H + v/>hT 

w w ‘^wjj [ 0*^0 o vT w wj 


‘"o S •''> ^ - "'w 'w • ^ 


<xi ) 
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If = constant, then the energy equation (v) can be 


written as 


+ v p H'T=^-{C(1-S)C(C»T-i-c/oST 

h O O Q W W wj GO WWW 


+ (1 - £) (pc)p T 


(xii ) 


On adding (xi) and (xii), we get 


K. + T - 

n 


{ 6 [< 1 -S„ 


€> O 


. s]} 

w w wj J 




?T f€r<'*-S >c T + S p c ill + (1-6.) (^C)„ ^ 
\ I w o o w w w JJ R 


This can be written as. 


2 

o^T T = T 


(xi i i > 


where 


£•{(1-3 )PC + sp c 1 + (1-6) (pc) 

woo W W W J 1 


O* = VpC+VpC_ 

o o o w w w 


Since i + j and v^- i + j 


theref ore 


where 


o* = O' i +0* J 

X y 


O' =rv O C + V pc 

X OX ^<y o wx ^w w 


O' =v pc + v pc 

y oy ^ ^ 
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In explicit form. Equation (xiii) is. 


dT <y & dT 

51 + (-y) ^ 

^t Sx iVj Sy 


P 

<— ) 


This equation is now written in the standard form suitable 
for the application of the OS algorithm as. 


^t 


+ 


U 


Sx 



(J2) 


(C> 


where 
U = 


O' 

X 




SCCi-S )p 
^ w o 


C + V pc 

o wx w w 

+ S P c y + (1-£)(pc)o 
w w w R 


and 

O' vpc+vpc 

y = — y = oy o vy w w 

e <:<1 - s ) p c + S ^ cJ + ( 1 - 8 )(|PC)„ 
w o o w w w *^8 

U and V are the composite velocities of oil and water along x 
and y directions respectively. 
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APPENDIX F 


(a) We show here that in a transformation from x-y coordinates 

to coordinates where ? is the streamline direction and r) 

is normal to it fSee Figure 4.19 (A)J, 

U T + V T s U' 

X y ? 

where 

U' = J 

At the origin of Figure 4.19(B), 

i+?yj and '^n = i + J 

^ -»v ^ 

= cos<jfr i + sin<^ j = — sin^ i + cos<^ j 

Using chain rule. 


y + T 

X ri 


r?. 


= T, ? + T 

? y » 


r). 


y < ■ y » y 

Substituting the expression in U T + V T , we get 

X y 


(U + V + (U + V Tjy ) 

= (U cos<j^ + V sin^) + (— U sin«^ + V cos4>) 


But cosi^ = 


U 


J U^+V^ 


, sirw^ = 


J U^+V^ 


Hence , 


U T + 

X 






“ Ju^ + T 


U* T„ 
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(b) Interpolation on a Two Dimensional Grid : 

The following method is used to find the temperature at 
an interior point in a two dimensional grid from the 
temperatures of the four surrounding nodes. It uses a 
linear two dimensional interpolating function. 

First (x*, y') is transformed to (ot, ft) fSee Figure 


4.19(E)> 




Q 


Ax 


- 1 


-1 i « 


Q 




1 


Now Tq 


2(y 


Q 




- 1 


Ay 


-1 ^ /Sq ^ 1 


4 

where is a weighting function. 

i=1 



The N.s are calculated as follows: 
1 


<1 - a > (1 - ft ) 

h 2- . Ng = 


(1 + a ^) (1 - ft ^) 


N 3 = 


(1 + Cq) (1 + /9q) 


' N 4 = 


(1 - o<q) (1 + /9q) 


4 


4 


APPENDIX G 

CALCULATING NODAL VELOCITIES FROM NODAL PRESSURES 
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From Darcy’s law CSee Section 4.7? Equation 4.6>, we get 


k k 


ro 


S p 

W 


Hence in component form, 
V = - d ^ , V 


i + 




Sy 


j ) 


OX 

Similarly, 


^x 
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Finite difference expression for the pressure 
the internal node points (2 i 2 j ^ 

follows s 


^Po 

dx 

l!£ 

dy 


p^ (i+1,j) — p^(i-1,j) 
2Ax 

p^ (i,j + 1) - Pjj<i»j-1) 
2Ay 


gradients 
M-1) is 


at 

as 


dp dp 

•^W ^W 


Similarly, expressions for 3fp" 

At the boundarieSf the expressions for pressure 


can be obtained- 
gradients are as 


follows , 


dx 


x=0 


dp. 

dx j x=X 



dp, 

dy|y=Y 


4 p(2,j) - p<3,j) - 3p(1,j) 

2Ax 

4 p(N-1,j) - p(N-2,j) - 3 p(N,j) 
2Ax 

4 p(i,2) - p(i ,3) - 3 p(i,1> 

2Ay 

m 

4 p(i,M-1) - p<i,M-2) - 3p<i,M) 
2Ay 



name of this program is NCW.F . IT SOLVES THE £-0 CONVECTION 
diffusion problem (UNIT VELOCITY ALONG XDIRECTION) USING 
OPERATOR SPLITTING TECHNIQUE . PARABOLIC PART IS SOLVED BY 
A.D.I. t ALTERNATING DIRECTION IMPLICIT ) TECHNIQUE . 

IT HAS THE FACILITY OF CHANGING THE B.C.S EASILY. 

dimension T1 (200.200 ) , T2(200.200) 

parameter (N0DESX»E6,N0DESY-£6) 

data T2/40000*0 . 0/ 

NHALF -13 


DO 40 I-1.NODESX 
DO 40 J-1 , NODESY 

IF( I .EQ.NODESX) T2(1,J)-0.0 
IFtJ.EQ.I) T2(I,J)-0.0 
IP(J.EQ. NODESY) T£tI,J)-1.0 
IFd.EQ.I .AND. J.GE. NHALF) T£(I,J)»1.0 
IFCI.EQ.I -AND. J.LT. NHALF) T2(I,J)-0.0 

) CONTINUE 

t^^^^^mv************ ************************************** ************ 

DO so I»1,N0DESX 
DO 50 J*l, NODESY 
0 TMI.J) ■ T2(I,Jr) 

XMAX»2. 0 
TMAX>*£ . 0 

PRINT * , 'WHAT IS PE , DELT? * 

READ * , PE, DELT 

DELX * XMAX /(NODESX-1 ) 

DELY * YMAX /(NODESY-1 ) 

TIME-O . 5 
DELT! - DELT 
DELTE ■ DELT 

NTIMESTEPS - IPIX (TIME^OELT) 

CALL PRINTDATA ( PE , T IHE . DELT , DELX . DELY , XMAX , YMAX . NODESX , 

1 NODESY, NTIMESTEPS) 

DO 1-1 .NTIMESTEPS 

CALL STEP1 (T1 , TE , NODESX , NODESY , DELT ) ,DELX) 

CALL STEPE t T 1 , TE , DELX, DELT, DELTE , PE , NODESX , NODESY , NHALF ) 

END DO 

CALL PRINT0UTPUT(T2, NODESX, NODESY ) 

END 

SUBROUTINE STEP1 ( T 1 , T2 , NX , NY , DELT . DELX ) 

DIMENSION T1 (200,200),T2(200,200) 

R-DELT/DELX 

DO I - 2 , NX-1 

DO J - E , NY-1 

T1(I,J) - R*TE ( I - 1 , J ) + 1 1 ■ 0-R ) -TE( I , J ) 

END DO 
END DO 

RETURN 

END 

SUBROUTINE STEPE (T1 . T2 , DELX . DELY . DELT2 , PE . NODESX , NODESY . NHALF ) 


dimension T 1 < 20 O,e 0 O>,TE(e 00 , 200 ) 

REAL INTERTEMP ( 200 ,200 ) , AR (200 ) . BR( 200) , CR(200 ) , FR( 200) , TEMPT ( 200 ) 

DELTADI - DELT2 / 2.0 

TEMPI « DELTADI / ( PE*DELX*DELX ) 

TEMPS - DELTADI / ( PE*DELY*DELY ) 

BLEFT ■ 1 + 2*TEMP1 
BRIGHT - -(20TEMP2) 

A - -TEMPI 
C - "TEMPI 
D - -TEHP2 
E * -TEMP2 

DO 1-1 ,200 
ARC I ) » A 
BR( I ) » BLEFT 
CRC I ) - C 
END 00 

DO J » 2, NODESY-1 

DO I« 2 , NODESX-1 

FR{ 1-1 )«{ 1+BRICHT)»T1 ( I, J)-D*T1 ( I, J+1 )-E*T1 (I, J-1 ) 

END DO 

FR( 1 ) • FR( I ) - A*T1 ( 1 ,.T) 

FR(N0DESX-2) - FR ( N0DES)C-2 ) - C»T1 ( NODESX , J ) 

CALL TRIDIA ( AR , BR , CR , FR , ( NODESX-2) , TEHP3 ) 

DO I«2, NODESX 

INTERTEMPd.J) ■ TEMP3(I-1) 

END DO 
END DO 

,,t********************************V¥m************m******m***m********m 

DO 20 1«1 , NODESX 
DO 20 J«1,NODESY 

1F( I .EQ. NODESX) INTERTEMP< 1 ..f)- 0.0 

IFtJ.EQ.I) INTERTEMPI I, J)«0.0 

IF( J.EQ.NODESY) INTERTEHP( I, J)-1 .0 

IFd.EQ.I .AND. J.CE.NHALF) INTERTEMP( I , J )«1 . 0 

IFd.EQ.l .AND. J.LT.NHALF) INTERTEMP d , J ) -0 . 0 

JO CONTINUE 

DO 30 I»1, NODESX 
DO 30 J«1,N0DE8Y 
30 TECIaJ) » INTERTEMPCI, J) 

BLEFT » 1 + S^TEMPa 

BRIGHT » -a^TEHPI 
DO 1-1,200 
AR( I ) « E 
BR(I) » BLEFT 
CRC I ) - D 
END DO 

DO I - 2 , NODESX-1 

DO J-2, NODESY-1 

FR( J-1 )-( 1 +BRICHT )* INTERTEMP (I, J)-C* INTERTEMP CI + l ,J)- 
» A*INTERTEMP<I-1 , J) 

END DO 

FR(1) - FR(1) - E*INTERTEHPd , 1 ) 

FR(NODESY-2) ■ FR(N0DESY-2) - D* INTERTEMP ( I , NODESY ) 

CALL TRIDIA ( AR , BR , CR , FR , ( NODESY-2 ) . TEMP3 ) 

DO J«2. NODESY-1 

T2d, J) - TEHP3(J-1 ) 

END DO 
END DO 


RETURN 

END 


i«pl ici t real ( k ) 
real infinity 

parameter (n«41 ,a»5, l*fi^ne«,biglim*0.1, itrnsli«*20) 
parameter (infinity»9 87 654381 .0, delta* 0.000001 ) 
parameter (nd*17) ! Domain dimensions, 

dimension po1(50,10),pwl(50>10) 
dimension po2( 50 , 1 0 ) , pwSC 50 , 10), po8oid( 50 , 1 0 ) 

dimension pw2old(50,10),sw(50,10),vdummy(50,10) 
dimension vox(50,10),voy(50,10),vwx{50,10), vwy ( 50, 1 0 ) 
dimension U{50,10),V(50,10) 
dimension told<50,10),tl(50,10),tm(50,10) 
dimension t2(50,10),t2old(50,10),si gma_t ( 50 , 10) 
dimension denoilC50,10),denwater<50,10) 
dimension A<61000),r(450),wm<450) 

integere4 Ike ep (2100), iw( 33 00), ivect (61000), j vect (61 00 0 ) 
inTeger»4 irn(61000),icn(122000),iglb(450,450) 


common /area1/delx,dely,deltl,Bi,ti 
common /area2/tbc_l , tbc_r 
common /area3/t1 , t m , t 2 , s i g»a_t , Kh 
Conversion factors . 


Conversion factor(p5i to pascal) 
Conversion f act or ( f ahrenhe i t to Celsius) 
Conversion f actorC megapascal to pascal) 
Conversion factoridarcy to meter-squared) 
Conversion factor(hour to sec) 


conv^psi = 6.894757e+3 
conv^f ah = 1 . 0/ 1 . 8 
conv^rapa = 1 . 0et6 
conv_darcy * 9.87e-t3 
conv^hour = 3600 

Parameter values Beta values are to be checked. 

2 eta_oil * 5 . 0e-6e ( 1 /conv_psi ) I unit == 1/pascal 
2 eta_water * 3 . 1 e-6* ( 1 /conv_ps i ) I 
beta_oil = 0 . 0004 1 ♦ ( 1 /conv_f ah ) 
beta_water = 0 . 0 0 04 1 ♦ ( 1 /conv_f ah ) 
eps i 1 on » 0 . 375 
k * 132 0*conv_darcy 


unit 

unit 

unit 

unit 

unit 


1 /pascal 
t /Celsius 
1 /Celsius 
dimensionless 




pw * 190.0 ♦ conv_psi / conv_mpa * unit * mpa 

tbc_l * 100.0 I unit » Celsius 

tbc_r = 50.0 1 unit = Celsius 

data tl ,t2, denoil, denwater/500*50.0,500m50.0,500<«900.0, 

1 500*1000.0/ 

do i*1 , n 
do j = 1 , m 

pwl(i,j) * pw 1 unit * mpa 

sw(i,j) » 0.2 1 unit * dimensionless 

call cap illarypres sure (sw(i,j),pc,ierror) 

po1(i,j) = ( pc/conv_mpa ) pw1(i,j) I unit = mpa 


pw1(1,j) * 260.0 * conv_psi / conv_mpa ! unit 
sw( 1 , j ) * 0 . 86 ! unit 

call capillary pressure (sw(1,j),pc,ierror) 
polClij) = C pc/conv_mpa ) + pwt(1,j) 1 unit 


mpa 

dimensionless 

mpa 


pw1(2,j) * 242.5 ♦ conv^psi / conv_mpa f unit » mpa 

s w ( 2 , J ) * .695 

call capillarypressure (sw(2, j), pc, i error) 

po1(2,j) = ( pc/conv_rapa ) + pw1(2,}) 1 unit = mpa 

pwtC3,J) = 225.0 * conv_psi / conv_mpa 1 unit * mpa 

sw(3,j) * .53 

call capillarypressure (sw(3, j>, pc, ierr or) 

po1C3,j) * ( pc/conv_mpa ) + pw1(3,j) 1 unit * mpa 

pwt(4,j) * 207.5 ♦ conv_psi / conv_mpa f unit = mpa 

sw ( 4 , j ) ■ . 365 

call capi 1 1 ary pressure ( swC 4 , j ) , pc , i error ) 

po1C4,j) = ( pc/conv_»pa ) + pwt(4,j) f unit =* mpa 


tHi.j) = I unit * Celsius 

* tbc_l I unit « Celsius 

t1(2,j) = ( t bc_l-t bc_r ) #3 . 0/4 . 0 + tbc_r 

t1(3ij) =* ( t bc_l +t bc__r ) / 2.0 » unit * Celsius 

t1(4,j) =* ( tbc_i-tbc_r ) *1 . 0/4 . 0 tbc_r 
♦ nd do 
end do 
do i*1 > n 
do j » 1 , » 

call oildensity ( po 1 ( i , j ) ♦conv_«pa , 1 1 ( i , j ) , deno ) 

call waterdensity (pw1(i,j) *conv_«pa , 1 1 ( i , j ) , denw ) 

denoil(iij) =* deno ‘ unit = k g/cub i c_»et er 

denwat er ( i , j ) = denw I unit » kg/cubic^neter 

end do 
end do 


delt * 0 . 0 1 ♦conv_hour 
tise * 1 . 0^conv_hour 
del X » 0.1 
del y » 0.8 
B i • 0.01 
ti = tbc_r 


I unit = sec 
I unit =» sec 
? unit « seter 
! unit » eeter 

! unit = Celsius 


sue » 0.0 
do 88 i * 1 , nd 


do 88 j*1 , e 
so * 1 . 0-sw ( i , j > 
sue - sue-t-so 

so^avgl = sue / float(ndee) 



ntieesteps * i f i x ( t i ee/de 1 t ) 

X eax * < n- 1 ) ♦de 1 X 
yeax * (e~1 )*dely 

call printdata(tieeidelt,deix,dely,xeax,yedx,n,a, 
ntieesteps) 

write(*, 1 1 15)nd,e,biglie, itrnslie,Bi , tbc_l , t bc_r 1 1 i 
foreat(' DOMAIN NODES ( X~D I RE CT I ON ) *',13/ 

' DOMAIN NODES ( Y-D I RE CT I ON ) *',13/ 

' BIGLIM *' ,F4 .8/ 

^ ITRNSLIM »' , 13/ 

' BIOT NO. (Bi) *',F6.3/ 

' INLET WATER TEMPERATURE *',F5.1/ 

'INITIAL FORMATION TEMPERATURE *',FS.t/ 

^ Ti «' ,F5. 1/) 


nh* 2^» 

c Giving values to IGLB. 
ni *0 

do 10 J * 1 , nh-t 1 

do 1 0 i* 1 , 1 “ j +1 ? lower 

n I » n 2 1 
irn(n 2 )*i+j -1 
icnCni >»i 

ICLB(irn(n 2 ),icn(n 2 )) = ni 
10 continue 

do 80 j« 8 ,nh -^'1 
do 8 0 i = 1 , l-j -^1 

n 2 *nx -^-1 ‘ upper 

irn ( nz ) «i 
icn( nz )*i^i ^1 

I GLB ( i rn ( nz ) , 1 cn < nz ) ) * nz 
20 continue 



lien = 2*"^ ^ 

lirn » ni + 5 

•typ* * ^ 
uf ■ ® ^ 
do 68 i-1 . 

jvect ( i )«icn( i ) 
cent inu« 


it*r 

V 

V 


. 0 

wsuK - 0.0 

osu« ® 0.0 
do 72 i-1 - " 

do 72 J*1 . ■ 

told( i . j ) = 

po2old( i . i ) 

pw2o 1 d ( i j ) ^ 
t2old( i . j ) 
po2 ( i j j ^ “ 

pw2 ( i < j ) 

t2( i . j ) 


T 1 ( i . j ) 

X po1(i<j^ 
* pw 1 ( i I ) ) 

t 1 ( i . i ) 
po 1 ( i » 3 ^ 
pw t ( i I 3 ^ 
t 1 ( i . 3 ) 


do 


2000 nt-l .ntimesteps 


delt 1 “delt 
n»ubt i«e5teps * 1 

i f ( nt . •q • ' ^ then 

delt 1 *delt/1 0 . 0 
nsubt imesteps * 1 0 

end i f 

ifCnt .eq. 2) then 

delt 1*delt/5 . 0 
nsubt i«esteps * 5 

endif 

if(nt .eq. 3) then 

delt 1 = delt/2 . 0 
nsubt iwest ep» * 2 

end i f 

do 1000 nt1-1 .nsubtisesteps 

bigeax ■ 0.0 

iterold * iter 


500 continue 


71 

c 


do 71 i*1 I ^ 
do 71 i*1 ^ » 

poSoldC i $ j ^ 
pw2old( i , 3 > 

t2old( i . 3 > 


po2 ( i » 3 ^ 
pw2 ( i . 3 ^ 
t2( i , 3 > 


Init ial ising . 

do 15 i=1»nz 

A(i)-0.0 

if ( i .gt . 1 ) go 

R ( i )*0 . 0 


ts continue 

c Bound.ry Conditions. 

do 110 if lag*' ' ^ 
if ( iflag.eq 1 > then 
ilower= 1 
iupper* 2*a 
istep*£ 
go to 108 

end i f 

if ( if lag . eq . 2 ) 


left b.c. 


? right b . c . 


ilow#r* S^Cn'-l 
iupp«r» 
i»t#p * 2 
go to t 09 
end i f 

if(ifl«g.eq.3) then 
ilower- 

iupper- ! botto* b e 

i»tep * t*m 
go to 109 
end i f 

if< iflag.eq.4) then 
i 1 ower* 1 

iupper* (n-t)*2*ii-1 ! top b.c. 

istep » 2^it 
go to 109 
end i f 

109 do 110 i»i lower , iupper, i«tep 
iree ®eod( i , 2) 

i f ( i re« . ne . 0 ) then 
in - < i+1 )/2 
else 

in * i/2 
end! f 

iy * mo6i in, «) 
if(iy.eq.O) iy * m 
ix « in/m + 1 

i f ( eod ( 1 n , • ) . eq . 0 ) ix « in/» 
r(i)* po1(ix,iy) 
r(i+1)= pw1<ix,iy) 
j 1 -i-nh 
} 2* i +nh 

if ( j 1 . It . 1 ) j 1*=t 

if(j2.gt.l) j2«l 

do 111 j* j 1 , j2 

A( IGLB< i , j ) ) * 0,0 

if ( i , eq . j ) A< IGLB( i , j ) ) * 1.0 

1 1 1 cont i nue 

i f ( j 1 . eq . 1 ) then 
j£-j2+1 
go to 112 
end i f 

if ( j2.eq. 1 ) then 
j1*j1+1 
go to 112 
endi f 
Jl-jl^1 
j2* j2^-l 

112 do 11 0 j* j 1 i j2 

A( ICLB( i'H , 3 > > * ^ 

ifd^l.eq. 3> A ( I CLB ( i i- 1 , j ) ) * 1.0 

1 1 0 cont inue 

do 200 i «2 , n- 1 
do 200 3*1, m 

c COMPUTING Ao, Bo, Co , Do,Eo, Fo AND Go VALUES 

cal 1 oi Irelpereeabi lity (sw(i, 3 ),kr 1 ,ierror) 
call oilviscosity (T2(i,J), viscoi 1 1 , i error ) 
ifCkrl . gt . delta) then 

terel * viscoill / C kr 1 ♦denoi I ( i , j ) ) 

else 

terel * infinity 
end i f 



C Computing th# values of thata towards th# #ast of tha noda . 
call oi Iralparaaabi 1 it y(sw( i-H . j ) , krE, iarror) 
call oilviscosity ( t E ( t , j ) , v isco i IE , i arror ) 
if(kr£ .gt. dalta) than 

taraE » viscoilE / ( krEadano i 1 ( i 1 , j ) ) 
alsa 

taraE » infinity 
andi f 

ifCtaral . It . ( infinity/E.O) .and.taraS. It . (infinity/E.O) ) 
tara - C t ara 1 i-t araE ) / E.O 
t aap -1.0/ tara 
a 1 sa 

t aap - 0.0 

and i f 

thata_a * ( k/apsi Ion ) ♦ (taap ) 

c Coaputing tha valuas of thata towards tha wast of tha noda . 
call oilralparaaability(sw( i-1 , j ) > krE, iarror) 
call oi 1 viscosity (tE(i-1 ,j), viscoilE, iarror) 
if(krS .gt. dalta) than 

taraE * viscoilE / ( krE^deno i 1 ( i*- 1 , j ) ) 
alsa 

taraS * infinity 
and i f 

if(t#ra1 . It . (infinity/E.O) .and. taraE. It , (infinity/E.O) ) 
tara » ( t ara 1 •i-t araS ) / E.O 
t aap » 1 . 0 / t ara 

alsa 

t aap » 0.0 
andif 

thata_w = ( k Zaps i 1 on ) ataap 

c Coaputing tha valuas of thata towards tha north of tha noda. 
if(j.aq.a) go to 401 

call oilralparaaability(sw(i,j+t),krE,iarror) 
call oilviscosity (tS(i,j-*'1), viscoilE, iarror) 
ifCkrE . gt . dalta) than 

taraE * viscoilE / ( krE^danoi 1 ( i , ) ) 

alsa 

taraE * infinity 
andif 

if(tara1 .It. (infinity /E.O). and. taraS.lt. (infinity /E.O)) 
tara » ( t eral +T: araE ) / E.O 
t aap *1.0/ tara 

alsa 

t aap * 0.0 
andif 

thata_n » ( k/apsi Ion ) ♦taap 

401 continua 

c Coaputing tha valuas of thata towards tha south of tha noda. 
if ( j . aq . 1 ) go to 40S 

call oi 1 r a 1 para alibi lity(sw(i, j-1 ),krS, iarror) 
call oilviscosity ( t E ( i , j -1 ) , v i sco i 1 £ , i arror ) 
if(krE .gt. daltu) than 

taraE * viscoilE / ( krS*danoi 1 ( i , j^l ) ) 

alsa 

taraE ■ infinity 
and i f 

if (ter* 1 . It . ( inf inity/a 0 ) end. ter*e . 1 1 . ( i nf i ni t y /a . (J ) ) 
ter* “ ( t er* 1 + 1 er*a ) / 3.0 
taap - 1.0 / tara 
alsa 

t aap - 0.0 
and i f 

thata_s * ( k/apsi Ion ) ♦taap 
40S continua 

if (j.aq. I) thata_s = thata_n 
if (j.aq.a) thata__n = thata_s 


than 


than 


t her 


t han 


c Co«puting the Ao values, 
anus * -theta_e 
denoM • deno i 1 ( i , j > #de 1 jcede I X 
Ao - anue/denoe 
c Coaputing the Bo values, 
so * 1 . 0-sw < i , j ) 
pc - po8( i , j )-'pwe{ i , j ) 

call sw_gradient ( pc«conv_apa, grad, terror ) 
teapi = (so*xeta_oi l-grad)/delt 1 
anua = thet a_e*f t het a_w 
denoa « deno i 1 ( i , j ) ♦de 1 x^de 1 x 
teap8 » anua/denoa 
anua * theta_n+thet a_s 
denoa « deno i 1 ( i , j ) ♦de 1 y ^de ly 
teap3 * anua/denoa 
Bo « teapti-teapEi-teapZ 
c Coaputing the Co values. 

Co » grad/deiti 
c Coaputing the Do values. 

Do = -theta_w/ (denoil(i,j )»delx^delx) 
c Coaputing the Eo values. 

Eo = ~t het a_n/ (denoil(i,j)*dely*dely) 
c Coaputing the Fo values. 

Fo = -t het a_s/ ( deno i 1 ( i , j ) ede ly^de 1 y ) 
c Coaputing the Go values. 

teap2 “ so*beta_oi !♦( 1 1 ( i , j )*told( i , j ) ) / deltl 
teap3 =■• ( tempi ♦po 1 ( i , j )+Co*pw1 ( i , j ) ) ♦ conv^apa 
Go » teapS ->• teap3 

c COMPUTING Aw,Bw, Cw,Dw,Ew,Fw AND Gw VALUES. 

call w«terrelperaeability<sw<iij),kr1,ierror) 
call weterviscosity (t2(iij)ivi sc waterl, terror) 
if(kr1 . gt . delta) then 

teral = viscwaterl / ( krl edenwat er < i , j ) ) 

else 

t era 1 * infinity 
end i f 

c Coaputing the values of gaaaa towards the east of the node, 
cal 1 waterrelperaeability<5w( i + t , j ) , krE, terror) 
call waterviscosity (t2(i + t,j) , viscwat er2, terror) 
if(kr2 . gt . delta) then 

teraE * viscwaterE / ( krE^denwat er ( t , j ) ) 
else 

teraE * tnftntty 
end t f 

1 f (term 1 . It . ( t nf init y/E. 0 ) . und . t eraS .lt.(tnfintty/E.O)) then 
term » < t er a t •♦•t eraE > / 2.0 
teap *1.0/ term 
else 

t eap * 0.0 
end! f 

gaaaa^e * ( k/eps i 1 on ) ♦ ( t eap 1 

c Computing the values of gaaaa towards the west of the node, 
call waterrelperaeability(sw( i“1 i j ) » krE, terror) 
call waterviscosity (tEii-t.j), viscwaterE, terror) 
if(krE gt. delta) then 

teraE = viscwaterE / ( k rE^denwat er ( i - 1 , j ) ) 
else 

teraE * infinity 
endif 

if (teral It . ( Inf ini t y/E. 0 ) . and . teraE . It . C inf init y/E. 0 ) ) then 
term » ( t era t •♦‘t eraE ) / 8.0 
t eap *1.0/ term 
else 


0 . 0 


T #*p ■ 

#nd i f 

ga«««_w * ( k/€p*i Ion ) ♦t 

C Computing th» values of gaaaa towards the north of the nod* 
if ( j . aq .») go to 403 

call wat#rralpar«aability(aw( i, j ^-1 ), krE, iarror) 
call watarviscoaity (tfiCi^ji^l), visewatarE, iarrorS 
if(krE .gt. delta) than 

taraE » visewatarE / ( krEadanwat ar ( i , j f t ) ) 
a 1 s a 

t araS « infinity 
end i f 

if ( taral .lt-(infinity/S.O).and . taraE It. (infinity /E.O)) 
tara * ( t era t -t-t araE ) / E.O 
t amp *1.0/ tara 
a 1 sa 

t aap ■ 0.0 
and i f 

gaaaa_n * ( k/apsi Ion ) ateap 
403 continue 


c Coaputing the values of gaaaa towards the south of the node, 
i f ( j . eq . 1 ) go to 4 04 

call waterralparaaability(sw(i,j-1),krE, ierror ) 
call water viscosity (t£(i,j-t) , visewaterE, i error) 
if(krS .gt. delta) than 

taraE * visewatarE / ( krEadanwat ar ( i , j«t ) ) 

else 

taraE * infinity 
endi f 

i f ( t arni . It . (infinity/S.O) .and. taraE .It. (infinity /E.O)) 
tara * ( t ar a 1 *tt ar aE ) / E.O 
t eap *1.0/ tara 
else 

t eap » 0.0 
endif 

gaaaa^s * ( k/epsi lon)ateap 
404 continue 

if (j.eq.1) gaaaa_s * gaaaa_n 
if (j.eq.a) gaaaa_n * gaaaa_s 


than 


than 


c Coaputing the Aw values, 
anua * -gaaaa^e 

danoa * denwat ar ( i , j )adal xadalx 
Aw * anua/denoa 
c Coaputing the Bw values. 

teapi * (sw(i,j ) a 2 ata_water-grad)/dalt 1 
anua * gaaaa_e+gaaaa_w 
danoa » danwater ( i , j )adal xadelx 
teapE • anua/danoa 
anua * gaaaa_n+gaaaa_s 
danoa » denwat er ( i , j ) •da 1 y»de 1 y 
teap3 * anua/danoa 
Bw » teapi '♦■t eapE-»'taap3 
c Coaputing the Cw values. 

Cw » grad/daltl 
c Coaputing the Ow values 

Dw * -gaaaa^w/ ( da nwat er ( i I j ) •de I X ada 1 X ) 

c Coaputing the Ew values 

Ew = -gaaaa_n/ ( dan water ( i , j ) •de 1 yade 1 y ) 
c Coaputing the Fw values. 

Fw * -gaaaa_s/ (danwaterC i , j )*delyadely) 
c Coaputing the Gw values. 

taapE = sw ( i * j ) abet a_wat era ( 1 1 ( i * j ) ** told(i,j)) / deltl 
taap3 » ( teapi apwt ( i , j )-^Cwapo! { i * j ) ) ♦ conv_apa 
Gw « teapS ^ teap3 



Ao * Ao*conv_«pa 
Bo = Bo^conv_»pa 
Co * Co*conv_iipa 
Do « Do^conv_iipa 
Eo » Eo«conv_»pa 
Fo * Fo*conv_«pa 

Aw = Aw»conv_«pa 
Bw * Bw*conv_«pa 
Cw * Cw*conv_apa 
Dw = Dw«conv_«pa 
Ew = Ew*conv_»pa 
Fw “ Fw^conv_apa 


if( j. 

aq . 1 ) 

t h*n 

Ao • 

0-0 


Bo * 

1 . 0 


Co * 

o 

o 


Do • 

o 

o 


Eo * 

-t . 0 


Fo * 

o 

o 


Co • 

0. 0 


Aw » 

o 

o 


Bw « 

1 . 0 


Cw • 

o 

o 


Dw • 

o 

o 


Ew • 

-1 . 0 


Fw » 

o 

o 


Gw * 
and i f 

o 

o 


if( j. 

aq . a ) 

then 

> 

o 

N 

0. 0 


Bo » 

1 . 0 


Co « 

0.0 


Do ■ 

o 

o 


Eo - 

o 

o 


Fo ■ 

-1 . 0 


Go ■ 

0. 0 


Aw * 

0.0 


Bw * 

1 - 0 


Cw » 

o 

o 


Dw » 

0 . 0 


Ew • 

o 

o 


Fw a 

-1 . 0 


Gw » 
and i f 

0. 0 



c A**igning Ao , Bo , Co , Do , Eo , Fo , Co and Aw , Bw , Cw , Dw , Ew , F w , Cw valua* 
c to tha ra*p#ctiv# A(i,j) slots and r(i>. 

nz1 » !K#aps track of tha noda in tha grid. 

nxa*£»nx1 ! Kaaps track of tha slot in A(i,j). 


A( ICLB(nza-1 . nxa-1 ) ) » Bo 

A( ICLB( nxa-l , nza ) ) * Co 

A( ICLB(nza~1 , nia-^-l ) ) - Eo 

A( ICLBCnza-^l , nxa«t +£♦») ) * Ao 

AC ICLBCnxa-l , nza-t ) * Do 

AC IGLBCnia-t ,nxa-3) ) * Fo 

r ( nza- t ) * Go 

A( IGLBC nza , nza ) > * Bw 

AC ICLBCnza, nza-1 ) ) « Cw 

A( ICLBCnia, nza-E) ) « Fw 



A ( I GLB( nz* , ) ) • Dw 
A ( I CLB( nz« , nz«-*-2 ) ) » Ew 

A ( I CLB ( nza , nz«4-2*« ) ) * Aw 
r ( nz« ) » Gw 


200 continue 

i f ( i t er . gt . 0 ) go to 55 

call ^aESa ( 1 , nz , A , 1 icn , i rn , 1 irn , i cn , uf , i k eep , i w , wa , if laga ) 
ifCiflaga.ne.O) print ♦i*iflaga-',iflaga 

i f ( i ter . eq . 0 ) go to 56 

55 continue 

call •a28b(l,nz,a,licn,ivect,jvect,icn,ikeep,iw,w«,iflagb) 
i f ( i f lagb . ne . 0 > print e^'iflagb *',iflagb 

56 continue 

call aaEScCliailicn.icn.lkeepfrfWaiatype) 

c Assigning solution r(i) to po2(i,j) and pw2(i,j). 
do 57 i=*1 , 1 
i rea *»od < i , 2 ) 
if ( irea . ne . 0 ) then 
in • ( ii-l )/2 
iflag • 1 
else 

in » i/2 
iflag * 2 
end i f 

i y » «od< in, a ) 
if(iy.eq.O) iy » a 
ix « in/a + 1 

i f ( aod( i n , a ) . eq . 0 ) ix * in/a 
i f ( i f lag . eq . 1 ) po2(ix,iy) * r(i) 
if ( iflag.eq.2) pw2(ix,iy) * r(i) 

57 continue 

c Updating sw(i,j) , denoil(i,j) and denwater ( i , j ) 
c (due to the pressure). 

do 58 i»2,n-1 
do 58 j«1 I a 

p ® po2(i,j) ♦ conv_apa 
t « t2( i , j ) 

call oildensity ( p, t , dens i t ynew > 
denoil(i,j) * densitynew 

p * pw2(i,j) conv_apa 

cal 1 waterdensity ( p. t , densitynew ) 

denwat er ( i , j ) * densitynew 

pc * po2 ( i 4 j ) '“pw2 ( i , j ) 

call watersaturat ion ( pc*conv_apa , sw( i . j ) , terror ) 

58 continue 

c Det era i nat i on of Darcy velocities, 
do 58 i* t , n 
do 59 j » 1 , a 

C 4 ll oilviscosity £ t 2 ( i , j ) , viscoi 1 , ierror ) 
call oi lr«lper«*«bi liTy(sw(i,j>,kro,i«rror) 
th*ti = ( k»kro ) /v iscoi 1 

call waterviscosit Y ( t2< i , j ) , viscw«ter. ierror ) 
call uat arralparaaabilityCawC i, j l.krw, iarror) 
gaaaa * ( k^krwl/viscwatar 

or i.eq.n or, j.eq.1 or j.eq.a) than 


if((i.#q.1 .or. i.oq.n) . ind . 
if ( i . #q . 1 ) ind#**! 
if ( i . #q . n ) ind#x*2 
c«ll •dgtqrad ( pograd , poS , i 
c«ll #dg#grad ( pwgrad , pwE, i 
ot#«px = conv_«pa ♦ pograd 
wtaapx = conv_«pa ♦ pwgrad 
( ( poac i , j + 1 )-po2( i 
( ( pw2( i , j + t )-pw2( i 


oteapy = 
wtaapy * 
go to 61 
andif 

i f ( ( j . aq . 1 . or . j . aq . a ) 

if ( j . aq . 1 ) index » 3 
if ( j . aq . • ) index * 4 
call adgagrad ( pograd , poa. i 
call adgagrad ( pwgrad , pwa , i 
otaapx = ( (po2( , j )-po2(i 

wtaapx a ( (pw2( i + 1 , j )-pw2( i 
otaapy * conv^apa ♦ pograd 
wtaapy = conv_apa ♦ pwgrad 
go to 61 
andif 

i f ( i . aq . 1 ) index * 1 
i f ( i . aq . n ) i ndax * 2 
call adgagrad ( pograd > po2 , i , j 
call adgagrad ( pwgrad , pwa , i , j 
otaapx a conv^apa ♦ pograd 

wtaapx a conv_apa ♦ pwgrad 

i f ( j . aq . 1 ) i ndax » 3 
i f ( j . aq . a ) i ndax » 4 
call adgagrad ( pograd , po2 , i , j 
call adgagrad ( pwgrad , pw2 > i , j 
otaapy = conv_apa ♦ pograd 

wtaapy * conv__apa ♦ pwgrad 

go to 61 
andif 

otaapx a ( ( po2 ( i + 1 , j ) -po2 ( i-l 
wtaapx = ( ( pw2 ( i + 1 , j ) -pw2 ( i-1 
otaapy » ( ( po2 ( i > j + 1 ) *po2 ( i , j 
wtaapy * ( ( pw2 ( i , j + 1 ) -‘pw2 ( i , j 


U-na.l .and, j.na.a)) than 


* j * n , a , index ) 
I j I n , a , index ) 


ij-1)) ♦ conv_apa)/ CEadaly) 
♦ conv_apa)/ (2adaly) 


and. ( i . na . 1 .and. i . na . n ) ) than 


# j / n , a , index ) 

# j > n > a , i ndax ) 

^ conv_apa)/ (2adalx) 
a conv_apa)/ (2adalx) 


.,n^a, index ) 
> n , a , i ndax ) 


^ n , a , i ndax ) 
,n,m, i ndax ) 


,j)) ♦ coriv_apa>/ {2*dalx) 
,j)) * conv_apa)/ (2tdalx) 
-1 ) ) ♦ coriv^apa)/ {2adaly) 
-1 ) ) » conv_apa)/ (2<dalY) 


cont i nua 


vox ( i , j ) 
voy ( i > j ) 
vwx ( i , j ) 
vwy ( i , j ) 

cont i nua 


that aaot aapx 
thata»ot aapy 
gaaaa^wt aapx 
gaaaa^wt aapy 


SOLUTION OF THE ENERGY EQUATION. 


CO *0 . 5a4 184.0 
cw « 4184 . 0 
Hr « 2412996 . 5 

Kh * 40. 0*4 . 1528a-3 


specific heat of the oil ( U--sac/Kg-“^ — -dtgC) 

specific haat of water (W-sac/Kg -dagC) 

heat “capaci t y / voluaa ( U-sac/ a* a 3 -dagC ) 
of the porous aadiua. 

Tharaal conductivity of the porous 
aadiua (U/a-dagC) 


do 6 0 i»1 > n 
do 60 j-1 , a 
so * 1 -sw C i , j ) 

taap! * apsilona ( so*co*danoi 1 ( i # j ) '•* sw( i , j ) acw^danwat ar C i ^ J ) I 
taap2 * ( I -apsi Ion) aHr 
sigaa^^t ( i , j I = teapl + taap2 

sigaa_x = vo x < i , j ) adano i 1 ( i , j ) ♦co vwx ( i , j I »danwat er I i , j } acw 



•g«« y “ voy ( i . j ) •denoll ( i , j )^C 0 + vwy ( i . J ) #r ( i , j)*cw 

\jii i) * slg«a_y / sig«a_t(i.i) 

^rroTOR ■ Te«P»ratur*_transfer is a subroutine that transr^rs the 
c predick , _,ture fro* the point q, which is upstreaa of t: K • point 
c tia**^*P This is the solution of the hyperbolic part 

^ to t ht po * ^ 


60 


if(i g« 2 i.le.n-1 .and. j.ge.£ and, j.le.a-n thien 

call stepi (U(i, j).V(i, j).i.i.ta(i, j)) 

*^!f((i.sq-1 O’- i.aq-«> »nd. i . ne . 1 and. i.ne.n ) th«.n 

call steplb ( U ( i , j > , i . 3 . t •< » > P > 

else 

tat i . 3 ) - tS( i , 3 ) 
end i f 
end if 
cent inue 


76 


O V 

, COmCTOR . Thi. p.r, .iU .oW. ... .-..He p.., 

C iaplicitly using A.D.I. and T.D.M. 

c*ll st^^p2 (n,*) 

, „pp..ln. Oil ..<■ -t.r — 
do 76 i=1 i n 
do 76 j*1 . ■ 

p = po2(i. j) * conv_apa 
t * t2{ i . i ) 

call oildensity (p, t , densitynew) 
denoiK i. i > - densitynew 

D « dw2( i , j ) * conv_apa 

call waterdensity ( p , t , dens it ynew ) 
denwaterti, j) - densitynew 
cont inue 

iter * iter + 1 
itrns = iter-iterold 
big * 0.0 
do 77 i^l.n 

do 77 j-1.a K f fnoRCi i)-po2old(i.3>> / 3 > > 

perron - 100. 0 « »bs ( ( po£ U . 3 » P . /pw2{i,j)> 

perrorS » ’ ^ 'i i-t2old( i , j 1 > / t2(i.j>> 

p.rror3 « 1 00 . Oaabs C ( t2( i . 3 i , 3^ ^ 

big * aaaxi (big.perrorl .perror ,p 

continue 

bigaax - ‘’!,^"*Lns le itrnslia) go to 500 

iftbig.gt .biglia .and. itrns. le.i 

do 40 i»1 . n 
do 40 3*^'® 
toldCi.j) = t1(i<3^ 
poi ( i , 3 ) * ^ ’ j ^ 

pwl ( i , 3 ) “ pw£< ' ' ^ ' 
ti ( i , j ) = t2( i , j ) 

cont i nue 


77 


40 


1000 


Co.pu..,.on of p.rc.o..,. poo-voU.. <PP»> 

sua * 0.0 
do 62 i»1 , nd 
do 62 j-1 . a 
BO » 1 . 0-Bw( i - 3 V 

sua ■ sua+BO 

so_avg * sua / float (nd*a) 


last 

Pf 


ppv = ^so_avgI - so_avg) ♦ 100 

c Computation of t i ■#~«vorag#d c#ntro~lino valociti#* for th# oil and 
c tha watar phasas. Oil valocity is calculatad at x»L and watar valocxty 
c is calculatad at x=L/2 (L * doaain_langt h y 
V_wsua » V_wsu« + vwx ( nd/a, (■'►t )/2) 

V_05U« a V_osua + vox ( nd, < a-t-l ) /a ) 

= V_wsu« / float (nt) 
o = V_osum / float(nt) 

i saa » mod ( nt , 5 ) 

i f ( i sea . aq . 0 ) than 

print *,'tiaastap «',nt 

print * , ' sub-t imast ap * ' , nt 1 

print ♦,'itrns'>itrns 

print 'big', big 

print a , ' sw ' 

do 107 )/a, )/a 

107 writa ( sw ( i , j ) , i » 1 , n ) 

1111 foraat(10<ax,f6.4)) 

print a , ' vox ' 
supfactor * 1 a-S 

do 108 j ■<«-«• 1 ) /a 1 ) /a 

108 writa (a/llia) ( ( vox ( i , j )/supf act or ) , i * 1 , nd , a ) 

Ilia format ( 1 0 ( ax , f 6 . 3 ) ) 

print ♦ , ' vwx ' 

supfactor » la-4 

do 1 08 j*(®+1 )/2, (•-►I )/a 

106 writa (a^lllB) !(vwx(i,j>/supfactor),i»1,nd,a) 

1113 foraat(10(ax,f63)) 
print ♦ , ' t ' 

do 113 )/a, i a^l )/a 

113 writa (a,1114) i 1 1 ( i , j ) , i * 1 , nd ) 

1114 format ( 10( 1 x,f7. 3) ) 

print ♦, ' ppv ' , ppv 
print ♦, ' ' 

andif 

£000 continua 

and 

subroutina stapt (U,V,i,j,tp) 
paramatar (pi*3.t415) 

common /araa1/dalx,daly,dalt1,Bi,ti 
common /area2/tbc_l , tbc_r 
Up » sqrt ( U^U+VaV ) 

if (U.ga.O .and. V.ga.O) indax * 1 

if (U.lt.O .and. V.ga.O) indax * Z 

if (U.lt.O .and. V.lt.O) indax » 3 

if (U.ga.O .and. V.lt.O) indax * 4 

i f ( abs ( U ) . 1 1 . 0 . 0 0 0 0 1 ) than 
tamp « pi/a . 0 
a Isa 

tamp » at an( abs ( V ) /abs( U > ) 
andif 

go to ( t 0 , £0 , 30, 40 ) , indax 
10 phi * tamp 

go to 50 

£0 phi = pi - tamp 

go to 50 

3 0 phi * pi tamp 

go to 50 

40 phi » a»pi - tamp 

50 continua 

daliata » Up ♦ daltl 

dalxq * dalzata a cos(piaphi) 

dalyq » dalzata ♦ sin(pi^phi) 

c Tampar atura_q calcuiatas tha tamparatura of tha point q which is ypstraam. 



c 


of thw point p i.».(i,j). 

call te«perature_q ( i . j . del xq, de I yq . t q, index ) 
t p » tq 

i f ( d# 1 z#t « . gt . de 1 X .or. delzeta . gt . dely ) then 
tp = ( tbc^l-ftbc^r )/2 . 0 
print ♦, 'ERROR IN STEP1' 
print ♦.i^j^delzeta.tp 
endif 
ret urn 
end 

eubroutine steplb (U,i,j,tp> 
coeeon /«r#a1/delx,dely,delt1,Bi,ti 
coeeon /area2/tbc_l, tbc_r 

coeeon /areal/ 1 1(50.10),te(50,10),t2(50,t0) ,sigea_t (50. 1 0 ) ,Kh 
delzeta * Uedeltl 
r = delzeta / delx 
if(U.ge.O.O) then 
t eep » 1 1 ( i-1 , j ) 
else 

t eep » 1 1 ( i-^1 , j ) 
endif 

tp * r^teep + C 1 -r ) et 1 ( i , j ) 
if(delzeta.gt.delx) then 
tp * ( tbc_l+tbc_r)/2 . 0 
print 'iRROR IN STEP1B' 
print ♦,i,j^delzeta,tp 
endif 
return 
end 

subroutine t eeperature^q (i,j,delxq,delyq,tq, index ) 

real n1 , n2 , n3 , n4 , net a 

coeeon /areal/ delx^delyjdeltt,Bi,ti 

coaaon /area3/t1 (50, 1 0 ) , T«( 50 , 1 0 ) . 1 2( 50 . 1 0 ) . sigaa_t ( SO . t 0 ) . Kh 
xp » (i-1)*delx 
y P * ( j - 1 ) ’■de 1 y 
xq ■ xp + delxcj 
yq a yp + delyc 
k i » 1 

i f ( index . eq . 1 .or. index. eq. 4) ki » 0 
i 1 » i- 1 +k i 
k j » 1 

if ( index . eq . 1 .or. index. eq. 2) kj = 0 
j1 « j“1+kj 

12 * i1+1 
j2 =. j1 

13 =* i 1 + 1 
j3 - jl+l 

14 » i 1 
j4 « jl+1 

x1 ■ (i1-1)*delx 
y1 » { j 1 -1 )*dely 
neta ■ ( 2* ( xq-x 1 ) /de 1 x ) - 1 
zeta “ { 2* { yq-y 1 ) /del y ) - 1 
n1 * ( 1 -net a ) » ( 1 -xet a ) /4 
n8 = ( 1 +net a ) • ( 1 -zet a ) /4 
n3 “ { 1 +net a ) • ( 1 ♦zet a ) /4 
n4 “ (1-neta)»( 1+zeta)/4 

tq - n1*Tl(i1,j1) ♦ n2*t1(i2,j2) + n3*t1(i3,j3) + n4*t1{i4,j4) 

return 

end 

subroutine step2 Cn,e) 
i ep 1 i ci t real ( k ) 

real interte«p( 50 .10),ar(50),br(50),cr<50),fr(50) ,Teap3<50 > 


r*« I arc ( 50 ) 

co««on /ar**! /del x , d*l y , d»lt 1 , Bi , t i 

co««on /area3/t 1(50.10),t«(50,10),t2(50, 1 0 ) , sig«a_t (50 , 1 0 ) , Kh 
daltadi ■ deltl/S.O 
cl ■ -2.0*dely*Bi 
c2 ■ -2 . 0*dely*Bi*t i 
call assign_bc ( ti» . n , ■ ) 

do j ■ 1 * ■ 

do i - 2 i n-1 

pe » 5 ig«a_t(i,j) / Kh 
teapi » deltadi / ( pe^del x ♦de I x ) 

t**p2 ■ deltadi / ( pe*del y ede I y ) 
a - -tempi 

bleft * 1 + B.OeteapI 
bright » 1 - 2.0*te«p2 

c » -tempi 
d » temp2 
e * temp2 
ar( i-1 ) = a 

br( i-1 ) * bleft 
cr< i-1 ) * c 

if ( j . It . m . and . j . gt . 1 ) then 
temp5 * tm( i , j+1 ) 
t emp4 * t a { i , j “ 1 ) 
endif 

i f ( j . eq - 1 ) then 
t emp5 * t m ( i , j 1 ) 

teap4 - tm(i,j+1) + cletm(i.j) - c2 

endif 

if ( j . eq . m ) then 

temp5 » t»(i,j-1) + c1*tm{i,j) - c2 

temp4 * tio(i/j~11 

endif . 

fr(i-l) * detempS + br ight *t m ( i , j ) e*temp4 

end do . 

(^) m fr(1) ~ ar(1)etm(1,3> 

f,.{n_ 2 ) - fr{n-2) - cr ( n-2 ) et m ( n, j > 

do i1 » 2,n-2 

arc( i 1 -1 ) • ar( i 1 ) 

end do ^ o -j \ 

call t r i d i a ( arc i br i cr # fr , n~ 2 # t e»pJ / 

do i* 1 # n-2 . X 

interteapC i***! f j ^ * teep3(i) 

end do 
end do 

call assign_bc( interte»p,n,*) 


1 


cJ o i * 2 I n t 

d o j * ^ ® 

pe = sigma_t ( i , j > / Kh 

tempi =« deltadi / ( pe*de I x ede I x ) 

Temp2 « deltadi / ( peedel y » de I y ) 

d * -Temp2 

bleft » 1 2.0etemp2 
bright - 1 - 2-0*temp1 

e » -Temp2 
a = tempi 


c " tempi , . , ■ 

s c» intertempC 1 + 1 , 3 

i nt er t emp < i , j ^ 
if ( 3 . ne . 1 • and . 3 . ne . m) 


> + a*iniertemp{ i-1 . j > 

then 


ar ( 3 -1 1 ® * 

br( 3 ) * bleft 
cr( j ) * <1 


end i f 

if ( j eq. 1 ) then 


+ 


bright* 



brd ) - 

cr(1 > “ ^ • 

fr(1) - > •*■ e*c2 

endif 

if (j.*P-*^ th»n 
br(«> ■ b lef T + d*c1 

ar(«-1) - d + • 

fr(*) *fr(ii) •*' d»c2 

e n d i f 
end do 

call tridi* («r,br, cr, fr, •, te»p3) 
do 3*1^* 

tad, j) - ttmplij) 
end do 
end do 

call asa i 9*^-^ ^ ^ ^ ^ ^ ^ 


return 

end 

subroutine ass ign^bc ( t , n # w ) 

di*tns ion t ( 5 0 ,1 0 ) 

coaion /areaS/tbc^l , tbc_r 

do 11 3*^ 
t(i,j) * tbc_l 
t(n,i) = tbc_r 
return 
end 

subroutine edgegrad ( grad , p , i , j , n index ) 
dimension p( 50 , 1 0 ) 

co««on /area1/delx,dely,delt1 >Bi >ti 
i f ( index . eq . 1 ) then 
pi * p( d j ^ 
p2 » p ( a / j ^ 
h a de 1 X 
go to 1 0 
end i f 

if (index. tq-2) then 
pi a p ( n, j 1 
pg a p( n-1 d ^ 
p3 * pCn-S/jl 
h a -delx 
endif 

if (index. eq. 3) then 
pi a p( i, 1 ) 
pg a p( i, 2) 

p 3 * p ( i I 3 ) 
h » dely 
endif 

if (index. eq. 4) then 
pi a p ( i , • 1 
pg a p( i , e-l ) 
p3 a p(i,«-2) 

h » -dtly 

endif ^ .v u ^ 

grad • (4 O^pS -“p3 --3.0»p1) / 

return 

grad « (pa-p1)/h 

return 

end 



